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ABSTRACT 


The free electron laser (FEL) is proposed to meet the Navy’s need for a speed- 
of-light high energy laser weapon capable of engaging a variety of targets including 
anti-ship cruise missiles, small boats, and theater ballistic missiles. A key attribute 
of FELs is good optical beam quality; in other words, they operate in only a few of 
the lowest-order transverse Gaussian modes. For weapons applications, a good mode 
quality is desired because it delivers the highest intensity on target ensuring a high 
level of lethality. A few higher-order modes can arise from the interaction of the 
electron beam with the optical beam, or from misalignments of the electron beam or 
resonator mirrors. High intensity on FEL optics can lead to mirror distortion due to 
heating and insufficient cooling of the mirror substrate. Mirror distortions, including 
astigmatism, can cause higher-order modes to appear affecting FEL performance. 
Therefore, it is important to quantify these higher-order modes because doing so 
uniquely identifies the optical field and may allow for corrective optics to single out 
the best modes for FEL lethality. 

This thesis will review free electron laser theory, and for the first time develop 
analytical solutions to quantify Hermite-Gaussian higher-order modes, develop a di¬ 
agnostic for modal analysis, and determine the tolerance limits on mirror distortions. 
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I. 


MOTIVATION 


Previous Naval Postgraduate School (NPS) classroom studies have shown that 
the PHALANX close-in weapon system, used onboard US Navy ships as a last tier- 
self-defense weapon, typically destroys an enemy missile at a distance of a few hundred 
meters from the ship. At that distance, NPS classroom studies have also shown a 
large probability that debris will cause collateral damage to the ship. A key factor 
in determining debris field collateral damage is time. An anti-ship cruise missile 
(ASCM) typically has a velocity of 1200 m/s. If it is detected at a nominal horizontal 
distance 1 of 15 km, then the ship has approximately 13 seconds to track and engage 
the missile. The probability of collateral damage stems from the fact that it takes 
several PHALANX rounds to destroy the missile. So, Naval interest in high energy 
lasers is clear: by engaging a target with a speed-of-light high energy laser weapon, the 
distance and time-to-engage increases considerably and the probability of collateral 
damage to the ship is reduced. 

In order to “kill” a missile, about one liter of missile material must be de¬ 
stroyed. This will cause sufficient structural damage resulting in missile breakup. 
The amount of energy required for material ablation is on the order of one megawatt 
over a surface area of 100 cm 2 for a few seconds. If an ASCAI, like the one described 
above, is detected at 15 km and engaged by a high energy laser, it would be destroyed 
in roughly five seconds at a distance of 9 km from the ship reducing the probability 
of collateral damage considerably giving the ship a large safety margin in time and 
distance. 

Several U.S. Navy laser programs are currently investigating the feasibility of 
employing high energy lasers on board ships. The free electron laser (FEL) possesses 

horizontal distance as a function of height of eye is given by D/nautical miles) = 
1.169,/Height of eye (ft) 
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many attributes that are in-line with the the Navy’s vision of an all-electric ship with 
no munitions: 

• High average power 

• Can be designed over a range of wavelengths (IR through X-rays) 

• Good optical beam quality 

• High operational reliability 

• No hazardous chemicals required for lasing medium 

• Virtually limitless magazine depth 

• Capable of engaging a variety of targets including ASCMs, small boats, and 
theater ballistic missiles 

A key attribute of FELs is good optical beam quality; in other words, they operate 
in only a few of the lowest-order Gaussian modes. A key attribute of FELs is good 
optical beam quality; in other words, they operate in only a few of the lowest-order 
transverse Gaussian modes. For weapons applications, a good mode quality is desired 
because it delivers the highest intensity on target ensuring a high level of lethality. A 
few higher-order modes can arise from the interaction of the electron beam with the 
optical beam, or from misalignments of the electron beam or resonator mirrors. High 
intensity on FEL optics can lead to mirror distortion due to heating and insufficient 
cooling of the mirror substrate. Mirror distortions, including astigmatism, can cause 
higher-order modes to appear affecting FEL performance. Therefore, it is important 
to quantify these higher-order modes because doing so uniquely identifies the optical 
field and may allow for corrective optics to single out the best modes for FEL lethality. 

This thesis will review free electron laser theory, and for the first time develop 
analytical solutions to quantify Hermite-Gaussian higher-order modes, develop a di¬ 
agnostic for modal analysis, and determine the tolerance limits on mirror distortions. 
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II. 


FREE ELECTRON LASER SYSTEM 


The free electron laser (FEL) system consists of many different components 
each of which involves very interesting physics. Figure 1 shows a basic systems di¬ 
agram of the operating FEL at the Thomas Jefferson National Accelerator Facility 
(.JLab). Let us walk through the operation of JLab’s FEL. A pulsed beam of electrons 
in a vacuum is produced using a photoinjector. The beam of electrons are accelerated 
to near the speed of light in the superconducting accelerator using radio-frequency 
(RF) energy. Highly relativistic electrons are then “wiggled” in the undulator pro¬ 
ducing laser light where it is amplified by extracting energy from the electron beam 
within the optical resonator. The electron beam’s remaining energy is recycled by 
having them enter the superconducting accelerator out of phase with the RF fields. 
Finally, the electron beam is directed to a beam dump. 

At the heart of the FEL system are relativistic “free” electrons, i.e., not bound 
to any atom or molecule. The electrons travel in a vacuum much like in a cathode ray 
tube. Conventional lasers operate at specific frequencies because electrons can make 
transitions only between discrete energy levels. The electrons of an FEL are forced 
to vibrate through a spatial periodic magnetic field. The frequency of vibration can 
be changed by adjusting the magnetic field or by changing the speed of the electrons 
which in turn changes the laser light frequency. 

In 1970 Dr. John Madey of Stanford University proposed what he termed 
a “free electron laser” using highly relativistic electrons to produce radiation in a 
magnetic undulator. Six years later Madev succeeded in demonstrating gain with a 
free electron laser with a 24 MeV electron beam and 5 m long undulator [2], 

As an illustrative example of an operating FEL, Table I shows the current 
parameters at JLab. At the time of this writing, JLab has the most powerful free 
electron laser in the world at 10 kW. The wavelength is in the infrared around 2-4 
//m. 
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Figure 1. Free electron laser systems diagram highlighting major components from 
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A. HOW DOES IT WORK? 

For an FEL to work we need [4]: 

• a beam of relativistic electrons, 

• co-propagating in an optical field, 

• through a series of spatially periodic magnets. 


The undulator establishes electron oscillations that are transverse to its overall longi¬ 
tudinal motion. In a reference frame moving along the FEL axis with the electrons, 
it would appear that they are simply oscillating and radiating in all directions. In the 
laboratory frame the radiation is confined to a narrow forward cone whose angle is a 
function of the electron’s speed, and is known as magnetic Bremsstrahlung radiation 
[2], The radiation emitted is highly Doppler shifted due to the relativistic nature of 
the electrons and the wavelength is given by 


A 


A, 


(I< 2 + 1 ) 

' 2 7 2 


(n. i) 
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Optical Resonator 

Output power 

10 kW 

Length 

32m 

Wavelength 

1. 1//111 

Waist radius 

0.6 mm 

Mode radius at mirror 

1 cm 

Quality factor 

25 (4% outcoupling) 

Rayleigh length 

0.5 m 

Max. intensity permissible on mirrors 

200 kW/cm 2 

Electron Beam 

Energy 

111 MeV 

Bunch charge 

135 pC 

Radius 

0.2 mm 

Peak current 

340 A 

Pulse length 

0.12 mm 

Undulator 

Period, A 0 

8 cm 

Length 

2.3 

Number of periods 

29 

Undulator parameter, K 

0.55 


Table I. JLab FEL Parameters from [3]. 


where A is the wavelength of the optical radiation, A 0 is the wavelength of the undu¬ 
lator magnetic field, 7 is the Lorentz factor, K is a dimensionless parameter that is 
proportional to the undulator magnetic field strength. We will derive this equation 
in the next chapter. 

The initial optical radiation given off by the oscillating electrons is known as 
spontaneous emission in the vernacular of conventional laser physics. The radiation 
from spontaneous emission is incoherent, meaning that, the wave packets emanating 
from each radiating electron are not in phase with adjacent wavepackets since each 
electron is randomly dispersed within the electron beam. Coherence evolves from this 
cacophony through mode competition with many passes through the undulator. The 
process of mode competition will be covered later. 
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At resonance, the electric field felt by the electrons oscillates at the same 
frequency and in the same plane as the electrons oscillating through the periodic 
undulator magnetic field. If the electric field and the electron’s transverse motion is 
parallel, electrons will lose energy and slow down since work is being done by them. 
If the electric field and the electron’s transverse motion is anti-parallel, electrons will 
gain energy and speed up since work is being done on them. The energy they lose or 
gain comes from the optical radiation field. The effect of electrons losing and gaining 
energy within each optical wavelength is a longitudinal bunching of electrons as some 
slow down and others speed up in each wavelength. The bunched electrons then emit 
radiation that is in phase and coherent. The intensity of the coherent radiation is 
much more powerful than the incoherent radiation [2], 

Two key performance metrics of the FEL can now be introduced: extrac¬ 
tion and gain. The fractional increase in optical radiation power in one pass of the 
undulator is defined as gain, 

P — P 

G = (II.2) 

O 

where the optical power, P, is measured at the end of the undulator and P 0 is the 
power at the beginning of the undulator. As the FEL converts electron kinetic energy 
to optical radiation energy, the fractional energy extracted from the electron beam in 
one pass of the undulator is defined as extraction, r/. 


Optical power 

V = —-- 

Initial electron beam power 


(11.3) 


The optical gain must be sufficient to overcome loss and outcoupling. A larger ex¬ 
traction means the FEL requires less average current to achieve high power. Typical 
values for extraction and gain are 2% and 30%, respectively. 


B. CONFIGURATIONS 

When the undulator is bounded within two reflective mirrors, the configura¬ 
tion is known as an optical oscillator or resonator. The optical radiation field that 
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couples to the electron’s motion can manifest itself through spontaneous emission. 
This radiation is collected within mirrors as shown in Figure 2. A partially reflective 
outcoupling mirror provides the optical feedback back into the undulator. Coherence 
develops through many passes of the optical field. The source of the optical radi- 



Figure 2. Free electron laser system in an oscillator configuration from [5]. 


ation field can also be an external laser in the case of an amplifier. An amplifier 
configuration makes use of a seed laser and typically a longer undulator. The longer 
undulator provides more gain and extraction in a single pass since only one pass is 
used. Coherence is established by the seed laser. Figure 3 shows a typical amplifier 
configuration. 


Electron Beam 


B field 


Accelerator ■ y 


Laser light 


c 


Seed Laser 


3 


Figure 3. Free electron laser system in an amplifier configuration from [5]. 
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C. COMPONENTS 
1. Source 
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Figure 4. Schematic of a typical photoinjector from [6]. 


An FEL requires high peak current in order for coherent radiation to grow in 
amplitude. The source for free electrons is a photoinjector. A drive laser, usually 
a solid-state laser, is pulsed at a photocathode target. Through the photoelectric 
effect, electrons are literally blasted from the surface of a cathode and accelerated in 
an electric field to the linear accelerator. Since the drive laser fires at regular intervals 
(on the order of nanoseconds), electrons travel in pulses. Electrons at the terminus 
of the photoinjector have energies from 5 to 10 million electron-volts (MeV). 

2. Accelerator 

The linear accelerator (see Figure 5) has two functions: (1) accelerate electrons 
from the photoinjector to its final velocity (energy ps 100 MeV), and (2) recover 
energy from recycled electrons. The electrons are accelerated as they travel through 
the cavities in phase with the radio-frequency (RF) fields. At JLab, electron bundles 
from the photoinjector are accelerated using a superconducting RF linear accelerator 
to minimize waste heat. Energy is recovered from recycled electrons and given back 
to the RF field since they reenter the accelerator out of phase with the RF field. 


8 












Figure 5. .JLab linear accelerator from [7] 


3. Undulator 

The undulator is a collection of spatially alternating permanent magnets that 
causes the electrons to wiggle transverse to their longitudinal motion. The transverse 
motion then couples to the optical electric field enabling work to be done on or by the 
optical radiation field. The magnets can be orientated linearly or helically as shown 
in Figure 6. Linear undulators produce linearly polarized light whose electric field 
oscillates in the same plane as the oscillating electrons. Helical undulators produce 
circularly polarized light whose electric field rotates along with the electrons enabling 
a coupling to occur. 

4. Beam Dump 

After energy is recovered from the recycled electrons in the linear accelerator, 
the electron beam is steered into a beam dump (see Figure 7). A beam dump is 
nothing more than a piece of notched metal that will absorb the remaining kinetic 
energy of the electron beam (on the order of several MeV). At at an average electron 
beam current of one ampere, about one megawatt must be absorbed. It is water- 
cooled on the non-vacuum side. 
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Figure 6. (a) Helical undulator resulting in circularly polarized light, (b) Linear 
undulator resulting in linearly polarized light from [8]. 

5. Optics 

As mentioned in Table I, the maximum permissible intensity on the resonator 
mirrors is on the order of 200 kW/cm 2 to prevent mirror damage to the mirror’s 
dielectric layers. In adequate or nonuniform mirror cooling may also warp the sapphire 
substrate thereby affecting FEL performance. FEL designs with a short Rayleigth 
length 1 have modes that expand so that the intensity on the mirrors is lessened. The 
entire optical system includes all the necessary beam directors to place laser energy 
on a target miles away from the ship. 


1 Ra,yleigh length is the distance the mode has to travel longitudinally from the waist so that the 

mode area doubles. Rayleigh length will be discussed in greater detail in Chapter IV. 
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Figure 7. Beam dump from [9] 
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III. 


FREE ELECTRON LASER THEORY 


We will now develop in detail the mathematics that describes the operation 
of the free electron laser. In order to produce laser light we must cause electrons to 
oscillate through a Lorentz force interaction, bunch the electrons spatially for coher¬ 
ence, and finally, have sufficient electron-to-optical mode coupling for amplification 
and gain. 


A. EQUATIONS OF MOTION 
1. Relativistic Lorentz Force 

In compact covariant notation, the relativistic Lorentz Force law 1 is 

e 


dp a dU a 
= m- 


= — F a ^U, 
c 


dt' dt' c~ “ /J (IIL1) 

where p a = ( 7771 c, p) is the four-momentum, 7 = 1 / \J 1 — v 2 /c 2 is the Lorentz factor, 
v is the electron’s velocity, m is the mass of the electron, £7 = ( 7 c, — 7 V) is the 
electron’s four-velocity, t' = t /7 is proper time, e is the magnitude of the electron’s 
charge, c is the speed of light, and F a ‘ 8 is the field-strength tensor with the components 
of the electric field, E, and magnetic field, B. Specifically, 

( 


jpa0 _ 




0 

~E X 

Ey 

~E Z 

E x 

0 

~B Z 

By 

Ey 

B z 

0 

-B x 

E z 

— By 

B x 

0 


(III.2) 


) 


Let us describe the electron’s trajectory in a helical undulator whose magnetic field 
along z is 

B u = B 0 (cos k 0 z, sin k 0 z, 0) (III.3) 

where B 0 is the magnetic field amplitude, k 0 = 27t/A 0 is the undulator wavenumber, 
and A 0 is the undulator wavelength, i.e., the distance for one complete cycle of the 


1 [10], Chapter 11 
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undulator magnetic field. The corresponding electric and magnetic field of the optical 
radiation field is: 

E r = E 0 (cos 0, — sin 0, 0) (III.4) 

B r = E 0 (sm 0, cos 0, 0) (III.5) 

where 0 = kz — cot + 0 , k is the optical wavenumber, uj = kc is the optical angular 
frequency, E 0 is the field amplitude in cgs units, and 0 is the optical phase. 

If we combine Equations (III.3), (III.4), and (III.5), Equation (III.2) becomes 

^ 0 — E 0 cos 0 E 0 sin 0 0 

E 0 cos 0 0 0 B 0 sin k 0 z + E 0 cos 0 

— E 0 simlj 0 0 — B 0 cos k Q z — E 0 sin 

0 —B 0 sin k 0 z — E 0 sin 0 B 0 cos k 0 z + E 0 sin 0 0 


j? a ) 3 = 




Using Equations (III. 1), (III. 6 ), and 


0 = - 
c 

7 = 




(III. 6 ) 

(m.7) 

(HI- 8 ) 


where (3 is dimensionless velocity and 7 is the Lorentz factor. We can write down the 
spatial components of the force from Equation (III. 1), 

^ = ~ — [Eg( 1 - !3 Z )(cos 0, - sin 0, 0) + d z B 0 (- sin k 0 z, cos k 0 z , 0)] 

at m 


d( 702 


me 

e 


g 

=- —[E 0 (/3 X cos 0 — fry sin 0) + B 0 (f3 x sin k 0 z — I3 y cos k 0 z)] (III.9) 

at. me 

where we have made use of f3± = f3 x i + f3 y j. Since (3 Z ^ 1 j (1 — f3 z ) —> 0, we then drop 
the terms proportional to E 0 ( 1 — ft z ) compared to l3 z B 0 . Integrating the transverse 
equation gives (assuming perfect injection 2 ) 

— 6 B 0 — 6 B 0 X Q 

/3 1 =--r-^-(cos k 0 z, sin k 0 z, 0) = —-——^(cos k„z, sin k„z, 0) (III. 10) 

7 me 1 k 0 zTT'ymc 2 


Perfection injection means that the electrons enter the undulator at precisely the correct angle 
such that its motion is exactly given by Equation (III.10) and constants of integration are zero. 


) 
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We can clean this equation up a little if we identify the dimensionless undulator 
parameter, K = eB 0 X 0 /('Imne 1 ) to write: 

—K 

/3 ± = -(cos k 0 z, sin k 0 z, 0) (III. 11) 

7 

We can now derive the energy exchange equation between the electrons and 
the optical field. With Equation (III.2), dt 1 = dt/j and p° = £/c = yrnc. from rela¬ 
tivistic kinematics 3 , we can write down the zeroth component, or energy component, 
of Equation (III. 1). For a = 0, we have 

^ = -'-F^U e 
dr c 


so that 


KAJ , 

7—17 me 
1 dV 1 


--(0,-E) • (7 c,-tv) 
c 


and 


d . . ey_, 

7—17 me) = -E • v 

'dtr' ’ c 

giving 

‘h = -J-0.E (IH.12) 

dt me 

Equation (III. 12) describes how the the electron exchanges energy with the optical 
field. The electron’s total energy, £ = ymc 1 . is a function of its rest mass and its 
speed (through 7 ). The undulator sets up electron motion in the transverse directions 
otherwise f3 ■ E would be zero. 

Substituting Equations (III.4) and (III.11) into (III.12), we arrive at 

eKE 

7 = t^ £ , oZ + T p) (HI. 1.3) 

7 me 


3 see [11], Chapter 1 


15 



Recall that 0 = kz — ut + (f>. If we define 

( = (k 0 + k)z — oot (III. 14) 

as the electron phase, we have 

eKE 

7= —SZtL cos(0 + 0) (III. 15) 

7 me 

We have arrived at an expression that describes energy exchange between the electron 
and optical field. Electron phase, (, is another way of describing the electron’s position 
along the undulator axis with respect to the co-propagating optical wave. Notice that 
for cos(( + 0) > 0 , the electron energy increases (absorption); for values less than 
zero, the electron loses energy (stimulated emission). 

2. Resonance Condition and Dimensionless Variables 
Introduced 

Since the electron phase, 0, is proportional to the electron position, z, let us 
define the dimensionless electron phase velocity as 

v=CL/c (III.16) 

where L is the length of the undulator. We can normalize the evolution of the 
electron’s motion to the length of the undulator by defining a dimensionless time, 
t = ct/L. It is convenient when possible to normalize our equations on scales that 
are relevant to the physics being explored. For example, longitudinal length can be 
normalized to the length of undulator. Thus, r takes on values from 0 —> 1 from the 
beginning of the undulator to the end of the undulator. In the following sections, we 
will normalize transverse directions, optical field amplitude, and current density. 
Returning to phase velocity, v can also be expressed as 

v=l (III. 17) 

where (.?.) = d(.. .)/dr is a derivative with respect to dimensionless time r. Phase 
velocity, zz. is a measure of how fast the electron is moving longitudinally within an 
undulator period. We can now derive the resonance condition. 
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The resonance condition results in optimum energy exchange between the elec¬ 
tron and the optical mode. This occurs when one wavelength of light passes over an 
electron as the electron travels through one undulator period. FEL resonance occurs 
when v = 0. Taking the time derivative of the electron phase gives us 

C = (ko + k)z — cot 

C = (k 0 + k)rhtc — co (III. 18) 


Substituting this expression into the definition of electron phase velocity gives us 


v = \(k 0 + k)/3*c — co\ — 

c 


(III.19) 


For v = 0, we have 


k = 


kofiz 


so that A = A 0 


(i - A) 

(i - A) 


A 


(III. 20) 


where we have used k = co/c. So now let us find an expression for /T in terms of 
physical parameters of the FEL. From Equation (III.11), /3j_ = K/j. and the definition 
of the Lorentz factor, we have 

1 

7 = 


7 2 = 


7 2 = 


1 

1 - PI ~ P s 


1 _ El _ Q2 

1 7 2 t J Z 


(III.21) 


Solving for the electron velocity gives us 


A = i / 1 - 77(^ 2 + 1 ) 


7 


(III.22) 
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For highly relativistic electrons employed in an FEL, 7 is typically 100, so we can 
make use of the binomial expansion 4 giving us 

A = \A- 7(^ + 1) 

ftz ~ 1-^(K 2 + 1) (111-23) 

Returning to Equation (III.20), we then have, 


A 


A, 


(.K 2 + 1 ) 
' 2 y 2 


(III. 24) 


Here is an expression for the wavelength of the optical mode based on physical pa¬ 
rameters that we can change (A 0 , K, 7 ): a clear example of the FEL’s tunability. It 
is now understandable why the wavelength of the light is so much smaller (by a factor 
of I/ 7 2 ) than the undulator period. This is physically plausible since the relativistic 
nature of the electrons causes them to see the undulator Lorentz-contracted by a 
factor of 7 and Doppler-shifted light from the electrons’ frame is then transformed to 
the lab frame reducing the wavelength by another factor of 7 . Electrons see resonant 
forces from undulator and the optical field. The end result is the optical field is at a 
much higher frequency then the undulator’s spatial frequency. 


3. Pendulum Equation 

We are now in a position to recast the equations in the previous sections into 
a form known as the pendulum equation. In this pendulum description, phase-space 
plotting becomes a natural and highly illustrative way of describing the electrons’ 
motion and interaction with the optical mode. Starting with the phase velocity, we 


4 The binomial expansion for small x is given by 


(1 + x) n = 1 + nx + n ^ 2! 1 ^ x~ + ■■■ 
see [12] pp. 14 for more details 
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have 


(III. 25) 
(III. 26) 


v = [(k 0 + k)/3 z - k]L =( 

v = (A : 0 + k ) fi z /. 

So let us now calculate an expression for ft... From Equation (III.23), 

A « 1-±(K* + 1) 

A = -)(- 2 )V 3 s( A ' 2 + 1 ) (m.27) 

Pz = - 3 (K 2 + l) (III. 28) 

/-yO 

But near resonance we have ( K 2 + l)/'/ 2 = 2A/A 0 . so Equation (III.28) becomes 

„ _ T(A' 2 + 1) 

Pz 9 

7 7 2 

0 ' 9 A 

= -v (III. 29) 

7 7o 

Returning to Equation (III.26), substituting Equation (III.29), and using A: = 27 t/A, 
we arrive at 

v = 

ZA = 

We can make the following substitutions and approximations: since the undulator 
period is much larger than the optical wavelength (i.e, A 0 >> A), and since the first 
term is smaller than the second by a factor A/A 0 , the first term is negligible. We define 
\ci\ = 47riVeA'|E , |L/7 2 ?7?.c 2 as the dimensionless optical field amplitude, L = N A 0 , the 
length of the undulator, dt. = Ldr/c , and substituting Equation (III. 15), we finally 
arrive at 

o 00 

v= \ci\ cos(C + (j)) =( (III.32) 

The pendulum equation describes the motion of the electrons in phase space, (Cw)- 


47tA 7 47T 7 

~^~ L+ Kj L 


(111.30) 

(111.31) 
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4. Phase Space 

Phase space diagrams are used in many physical problems. For example, in a 
simple pendulum, the angular displacement, 9(t), and its angular velocity, v(t), can 
be thought of coordinates of a two-dimensional phase space. As the motion of the 
pendulum evolves, the point given by ( 6 , v) will move along a path in phase space. 
For different initial conditions, the motion will be described by different paths. Phase 
space plots can be divided into two categories: closed orbit or open orbit paths. 

For the pendulum, closed orbit paths result when the pendulum is bound 
about the stable fixed point. Open orbit paths occur when there is a large enough 
angular velocity to cause the pendulum to go over the top. 


Pendulum Phase Space Paths 



0 


Figure 8. Various phase-space paths for a simple pendulum. The thin green curves 
are examples of closed orbits. The dotted red curves are examples of open orbits. 
The separatrix is the heavy dark curve. 


Figure 8 illustrates the phase space paths of a simple pendulum for 10 different 
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initial conditions. The phase space paths are given by 5 

v 2 = v 2 + 2(cos 9 — cos 9 0 ) (III.33) 

The separatrix is found when 9 0 = 7r and v 0 = 0 in Equation (III.33) giving us 

v 2 = 2(cos 0 + 1) (III.34) 

The separatrix separates locally bounded motion from locally unbounded motion and 
passes through the unstable fixed point (9 0 = 7r, v 0 = 0). 

Equation (III.32) is equivalent to the simple pendulum if |a| and 0 are constant, 
as is approximately the case for low gain in the FEL. The electron phase, £, can be 
viewed as the electron’s position within an undulator period while a change in phase 
velocity is proportional to the electron’s energy (recall v oc 7 / 7 ). The dimensionless 
optical field amplitude, |a|, is a measure of the field strength. A larger |n| causes 
the phase space paths to evolve much faster (analogous to a larger \fgjl for the real 
pendulum). A larger |a| also causes the separatrix to increase in size (along v). The 
FEL separatrix is given by [13]: 

v 2 = 2|a| [1 + sin(£ T 0)] (III.35) 

Figure 9 shows the evolution of 1000 electrons in phase space all with an initial 
v 0 = 0 (on resonance) and \a\ = 1. Notice that half the electrons gained energy from 
the optical field (increasing values of v) and half the electrons lost energy to the 
optical field (decreasing values of v). Notice how the electrons are becoming bunched 
along £. Since there are an equal number of electrons gaining and losing energy, the 
overall gain is zero. Gain is the fractional change in the power of the optical field 

“Dimensionless pendulum equation given by 0 = — sin 6. Multiply both sides by v to get v v= 
—v sin 8, which is equivalent to d[v 2 / 2 — cos 6\/dt = 0, giving us Equation (III.33). 
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Figure 9. Phase space evolution for 1000 electrons with all with initial v = 0. 



Gain 


inniiiiinmiiiimiiiiiii 


over a pass through the undulator. In symbols, we write: 


G 

and 

P 


P-Po 

P 0 


(III.36) 


(III.37) 


where the power, P is measured at the end of the undulator and P 0 is the power at 
the beginning of the undulator. We will discuss gain in detail in the next section. 

In Figure 10, we start out with an initial u 0 = 1.5 and |a| = 5.0. There is 
an imbalance in the number of electrons that are gaining and losing energy from the 
optical field. In this particular case, there is a clear bunching of electrons near ( = it 
resulting in gain. 


B. ELECTRON BEAM AND OPTICAL FIELD INTER¬ 
ACTION 

In the previous section, we derived the FEL pendulum equation which tells us 
the phase space evolution of electrons in an undulator period. Phase space evolutions 
that ended with good bunching at the appropriate phase resulted in gain (see Figure 
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Figure 10. Phase space evolution for 1000 electrons all with initial v = 1.5. 



10). The upper picture of Figure 11 illustrates bunching in an optical wavelength for 
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Figure 11. In the upper picture, electrons have an initial v 0 = 0 which results in 
bunching but no gain. In the lower picture, electrons have a an initial v 0 > 0 resulting 
in good bunching and positive gain from [4], 


u 0 = 0 or on-resonance case. At time the electrons start out evenly spaced within 
a section of the electron beam one wavelength of light long. At time f 2 , the electrons’ 
positions have evolved so that the number of electrons which give energy to the optical 
field (gain) equals the number of electrons which take energy from the optical field 
(absorption). The result is zero net gain but good bunching. In the lower picture, at 
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time t\ evenly spaced electrons have an initial non-zero phase velocity, i.e. v > 0. At 
time t‘> the electrons’ positions have evolved so that there are more electrons giving 
energy to the optical field than absorbing it. The result is positive gain, i.e., the 
optical field grows. 

1. Maxwell’s Equations 

Let us quantitatively describe the interaction between the electron beam and 
the optical field. We begin with Maxwell’s Equations in cgs units 6 7 : 


V-E 

= 47T/9 

(III.38) 

VB 

= 0 

1 <9B 

(III.39) 

V x E 

c dt 

4tt IdE 

(III.40) 

V x B 

= —^ ^ - TT 

c c at. 

(III.41) 


where p is the charge density, E is the electric field, B is the magnetic field, and J 
is the current density. Equation (III.39) automatically implies that B is the curl of 
some other vector, A, the magnetic vector potential'. We then have 


VB = 0=4>B = V xA 


(III.42) 


If we substitute Equation (III.42) into Equation (III.40) we get 

1 d 

V x E = — — (V x A) 
c at 

i d 

V X E + - — (V x A) = 0 
c at. 

„ / idA\ 

Vx ( E + t ~m) “° 

Recall from elementary vector calculus, if the curl of some vector is zero, then that 
vector is equal to minus the gradiant of a scalar 8 , <f>, the electric potential. We can 


6 see [10], Chapters 6 and 11 

7 [14], Chapter 5 

8 Ibid 
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then write 


1 d A \ 1 d A 

V x ( E + = 0 => E + = -V$ 

cat) c at 


Therefore, 

1 d A 

E = -V$-— (III.43) 

c dt v ' 

Equations (III.42) and (III.43) describe the electromagnetic radiation in terms of their 

potentials, A and If we substitute them into Equation (III.41) we arrive at 9 


VxVxA=-J|~[ -VT - ] 

c cat 


V(V • A) — V 2 A 

a 1 9 2 A 
V A - 


47T 


dt 2 


-J --V(— 

: c at 

„ . ia$ 

V-A+- — 

c at 


iaA\ 

c dt J 

i d 2 A 

: 2 dt 2 
47T , 


(III.44) 


For completeness, if we substitute Equation (III.43) into Equation (III.38) we get 


i dA 

-V$-— 

c dt 

1 d 


= 47 Tp 


V 2 T + - —(V • A) = -4t rp 
c dt 


(III.45) 


Equations (III.44) and (III.45) fully describes the optical field. By a judicious choice 
of Coulomb gauge, V • A = 0, Equation (III.44) becomes the full wave equation in 
terms of transverse current 10 

1 <9 2 A 47r , 


V A 


(III.46) 


e 2 dt 2 c 

We will use this equation to determine the electron beam to optical mode interaction. 

2. FEL Wave Equation 

For our helical undulator given by Equation (III.3), a solution to 

1 d 2 A 4t r , 


V A 


dt 2 




(III.47) 


9 Vector Identity 8, [14] 

10 [10], Chapter 6 
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can be written as 


A(x, t) 




(III.48) 


where E(-x,t) is the complex optical amplitude, a = kz — cut is the carrier wave, 
i = (—*,1,0) is optical field’s polarization vector for circularly polarized light, and 
x = xi + yj + zk is the position vector. Substituting Equation (III.48) into (III.47), 
we get the following 


d 2 A _ ee ia d 2 E 
dx 2 k dx 2 

d 2 A ee ia d 2 E 


dy 2 k dy 2 


d 2 a 

ee m 

'd 2 E 

dE 


dz 2 

k 

_ dz 2 

+ 2/A:—— 
dz 

- k 2 E 

d 2 A 

ee ia 

'd 2 E 

dE 


dt 2 

k 

dt 2 

1 

3 

' <S> 

CM 

1 

- uj 2 E 


(111.49) 

(111.50) 

(111.51) 

(111.52) 


For a laser beam, we can make a few simplifications. If we assume that optical wave 
amplitude and phase are slowly varying along the ,2 axis, second derivatives in z and t. 
can be dropped compared to first derivatives. Substituting these back into Equation 
(III.47) and using uj = kc , we arrive at 


ee 


A: 


o 0 ., I d Id 
Vj_ + 2?,A: ( ——I— — 


E = -^ J 

c 


fz-cuy- -- (IIL53) 

where = d 2 /dx 2 + d 2 /dy 2 is the transverse Laplacian operator. We can simplify 
the expression further by a change of coordinates, u = z — ct , which follows the light 11 , 
and recalling dimensionless time, r = ct/L , then 


d_ 

dz 

d 


dt 

, d Id 

so that, ——I-— 

dz cdt 


du d dr d d 

dz du dz dr du 

dr d du d c d d 

dt dr dt du L dr ° du 

d led d 

du c L dr du 

\_d_ 

L dr 


(III.54) 


11 


Coordinate transformation is also known as method of characteristics. See 


[15] 
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We have effectively collapsed z and t into one variable r. Multiplying Equation 
(III.53) by ke*e m , making use of Equation (III.54), and e* • e = 2, we now have 

V* + 2 ik (4^-1 E = - 27 T-er ia J ± • e (III.55) 

\L or J \ c 

We can express the electron current density as a sum over all electrons 

Jj_ = —ec E /3±S^(k — x;(t)) (III. 56) 

i 

where e is the magnitude of electron’s charge, x; is the electron’s position at time t, 
and S ^ is the three-dimensional delta-function. Equation (III.11), /3j_, can also be 
expressed as 

= »{- — iee- ik ° z } (III.57) 

7 

where ti signifies taking the real part of the complex argument. Combining these 
equations, Equation (III.55) becomes 

+ 2 ik ( — —— I E{ r, r) = -4tt ikeKp < -— > (III.58) 

\ L ot J \ 7 

where r = xi + yj. We have made use of Equation (III.14) and 

p(x,t)= / d)^(x — x; (t))dV (III.59) 

Jv 

The average < ... > represents an average over sample electrons in a volume element. 
If we multiply Equation (III.58) by —AniNeKL 2 /j 2 me 2 k, recall the dimensionless 
optical field amplitude, \ci\ = 4.irNeK\E\L/^rnc 2 , and define a dimensionless current 
density j = 8 x '.\7 2 1\ 2 1//>/- ‘‘me 2 , we have 

^ a(r,T) =< j( * > (III.60) 

when 7 & y 0 , a reasonable assumption for all the electrons during the length of the 
interaction. Looking at the factors multiplying the transverse Laplacian suggests we 
define dimensionless transverse coordinates (where the tildes refer to the dimensionless 
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variables): 


(111.61) 

(111.62) 

We can now recast the Laplacian operator into a dimensionless Laplacian operator, 
V^_, since 

8 2 d 2 _ k Id 2 d 2 

dx 2 dy 2 2 L dx 2 dy 2 

The FEL wave equation now becomes, 

+ "(f : r) < :i< K > (III.64) 

where a(r,r) = |a|e®^, the complex dimensionless optical field amplitude, measures 
the optical field strength in the interaction region; the dimensionless current j mea¬ 
sures the interaction between the electron beam and the optical mode; and < e ^ > 
measures the amount of electron bunching. If j is large ( j » 7r), the optical field a 
changes rapidly so we have high gain. If j is small ( j « i r), we have low gain. 




3. Gain 

Thus far in our development of FEL physics, we have studied free electrons 
interacting with laser light through a Lorentz force interaction in the undulator. We 
have described, through the pendulum equation, the process of physically bunching 
the electrons in phase space leading to coherent radiation. Finally, we have coupled 
this radiation, through Maxwell’s equations, with the electron beam co-propagating 
in the undulator with the optical beam (see Equation (III.64)). 

We must now discuss the modification of the field by FEL amplification. In 
weak fields, |«| << 7r and low gain j « n, we can apply a perturbation expansion of 
the pendulum equation resulting in an expression for optical mode amplication [13] 
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where the filling factor, F, is the ratio of areas of the electron beam to the optical 
mode. Gain is proportional to the dimensionless current density, j. At r = 1, the 
final gain is 




2 — 2 cos(^ 0 ) — v 0 sin(i/ 0 ) 


(III.66) 


Equation (III.66) is shown in Figure 12. Notice the anti-symmetric nature in v 0 


Gain vs. v ati = 1 

0 



Figure 12. Gain curve plotted at r = 1. 


and zero gain at v 0 = 0. The gain peaks at v 0 ~ 2.6 with peak absorption at 
v 0 ~ —2.6. When the electrons are at phase velocities slightly above resonance, 
the amount of energy transferred from the electron beam to the optical mode is 
maximized meaning that the optical mode is amplified by the electron beam. In 
an operating FEL, electrons that initially go through the undulator spontaneously 
emit in a spectrum (Figure 13a) centered around the resonant frequency. The weak 
field low gain curve (Figure 13b) narrows the spectrum over many passes (Figure 
13c). Over many passes, the laser spectrum narrows due to mode competition. The 
top of Figure 13d shows a narrow linewidth after 2500 passes requiring only a few 
microseconds of FEL operation. 
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-6 V 6 

Figure 13. FEL coherence as it evolves from a broadband spectrum, (a) Broad¬ 
band spectrum from spontaneous emission, (b) Low field gain curve, (c) Narrowing 
broadband spectrum over many passes, (d) Final power spectrum showing narrow 
linewidth. 
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IV. 


OPTICAL THEORY 


In the previous section, we described how oscillating electrons co-propagating 
with the optical field gave rise to amplified laser light. We will now restrict our 
discussion to the the laser light itself as it propagates from the interaction region in free 
space. We are interested in the transverse optical mode structure because the burn 
pattern (ex \a(x, y) | 2 ) will determine the intensity delivered to a beam optic or a target. 
Free electron lasers (FEL) primarily operate in the fundamental Gaussian mode, but 
a few higher-order modes are also possible. It has been shown in simulations and 
experiments that a few higher order modes can develop when the FEL is operated 
with mirror shifts and tilts, electron beam shifts or mirror asymmetry. Building on 
previous work [5], we will develop the mathematics for Hermite-Gaussian modes and 
a means of finding the coefficients for a given waveform. The modal analysis described 
here will be applied to FEL mirror distortions which is being investigated for the first 
time. 


A. OPTICAL WAVE EQUATION 


We first begin our discussion of optical theory by describing light propagation 
when no sources are present. Recall from the previous chapter, the FEL wave equation 
in the undulator is 


i ~ 2 d 

-tV + TT" 


«(r,r) =< — je > 


4 ’ dr_ 

when no sources are present, i.e., j = 0 and Equation (IV.1), reduces to 


(IV.l) 



d ' 

+ 4 * 77 - 
OT 


ci( r, t) = 0 


(17-2) 


where a = |a|e'V, the complex dimensionless field amplitude, measures the optical 
field strength. In terms of the FEL system, the following description is applicable 
outside of the undulator or when the electron current is zero inside the undulator. 
In our definitions of dimensionless variables x, y, and r, L represented the length of 
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the undulator. It will now represent the length that the optical mode travels. Our 
task now will be to find solutions of Equation (IV.2) to investigate transverse mode 
properties. 

B. FUNDAMENTAL MODE 

A form for the solution of Equation (IV.2) in lowest order can be written down 
as 

ci(r,r) = a Q e (IV.3) 

where a( r, r) is the dimensionless complex field amplitude proportional to the electric 
field, f 2 = x 2 + y 2 , and a 0 is the real dimensionless optical amplitude at x = Q,y = 
0,t = 0. Functions p(r) and q(r) are yet to be determined. Rayleigh length, z 0 = 
Z 0 /L , is the longitudinal distance from the mode waist at which the optical area 
doubles, i.e., where the mode radius exceeds the value at the waist by a factor of y/2. 
It is a measure of how the optical mode expands (see Figure 14). Equation (IV.3) 


a 


b 


0^1 

Figure 14. Cross section of a Gaussian mode propagating from r = 0 —> 1. The mode 
with a smaller Rayleigh length expands more quickly, (a) z 0 = 0.20, t w = 0.0 (b) 
z 0 = 0 . 8 , r w = 0.0 
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is the dimensionless analogue of a form used by [16], et al. Since we are solving a 
second-order differential equation, we must specify initial conditions to determine a 
unique solution. We stipulate that at r = 0, a( r, 0) = ci 0 e s °, an initially transverse 
Gaussian wavefront. The initial conditions then require: 

p( 0) = o 

q( 0) = 1 (IV.4) 

Our job will be to find p(r) and q(r). 

We now take derivatives of Equation (IV.3) resulting in: 


d 2 < 


-2 4£ 2 


= U- + 


dx 2 \z 0 q z 2 q 2 


(IV.5) 


d 2 a _ f— 2 + 4 y 


dy 2 


z 0 q z 2 q 2 


(IV.6) 


da ( ° x 2 Q y 2 q\ 

d^~ \ P+ J^ + J^) a ( ' } 

(Recall (.?.) = d(...)/dr.) Combining Equations (IV.5), (IV.6), (IV.7) into (IV.2) 
and dividing through by a, we are arrive at 


2 , 4 ^ 2 2 , # 2 , ,. / ° , x 2 q y 2 q _ 

— — -h Tyy — --h 9 + 4% I — P + -—- + -—- | — 0 

z a q z 2 q 2 z 0 q z 2 q 2 \ z 0 q 2 z 0 q 2 


Making use of f 2 = x 2 + y 2 , we have 

.2 / 1 


z 0 q \z 0 


— + / q I = \ — + ip 


z 0 q 


(IV.8) 


(IV.9) 


For this equation to hold for any r 2 at some particular time r, both terms in paren¬ 
theses must be identically zero. Solving the left hand side for q(r) 


q = — 


i 

Z 0 


q(r) = 


-dr' 

Z n 


IT 

= q 0 + — 

Zn 


(IV. 10) 
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Applying the initial conditions, Equation (IV.4), we find 


, \ lT 

q( T ) = i + — 


(IV. 11) 


Solving the right hand side of Equation (IV.9) for p(r) and substituting Equation 
(IV. 11), we find P(r) to be 

(IV. 12) 


p = lZn 


+ 


T 


i 2 + r 2 z 2 + t 2 

Integrating with respect to r and applying initial conditions, we have 


p{T) = J tw^ + s^ )dT ' 


(IV. 13) 


n- I n- 2 _L 2 ? 

p(r) = i arctan — + In ' 0 


(IV.14) 


Substituting Equations (IV.11) and (IV. 14) back into Equation (IV.3), we find 


a( r, t) = a 0 e 


-i arctan In 


(IV.15) 


We can clean up this mess if we make a few substitutions. From Equation (IV.11), 
we have 


1 

(l(r) 


1 + V 

Zo 


T 


— I- 


Z 0 ~\~ — 

0 Zo 


Z 0 -j- — 

0 Zo 


With the benefit of hindsight, we can make the following definition: 


w 2 (t) = z 0 +^~ 



(IV.16) 


(IV. 17) 


Equation (IV. 17) represents the optical mode radius as shown in Figure 15. The term 
z 0 is the dimensionless analogue of the dimesional Rayleigh length. The derivation 
has assumed that the beam waist occurs at r = 0, that is, (7(0) = w 0 = But 
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the beam waist may be positioned anywhere along r. It is a trivial matter to make 
a time translation in the equations such that, r —> (r — r w ). With Equation (IV.17), 
we can re-express l/q in terms of w and w 0 . Equation (IV.16) now becomes 



(IV.18) 


Figure 15. Cross section of a Gaussian mode as it propagates from r = 0 —> 1 
illustrating the optical mode radius. r w = 0.3, z 0 = 0.3 


Returning to Equation (IV. 15) and substituting Equation (IV. 18), the optical 


field becomes 1 


w n 


=2 . * 2 . 
I- 


a(r,T) = a 0 , . . e * 2 (^)e 


-i arctan 4- 


W[T ) 


(IV.19) 


Finally, we can write down the fundamental mode of the optical field as 


Amplitude 


w n 


a(r,T) = a 0 — . e * m e 


H-r) 


-arctan — 

Zo 


WIT) 


Phase 


(IV.20) 


1 The square root term can be expressed more compactly as w 0 /w 
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where the mode waist is at r w = 0. For any arbitrary location of the mode waist, we 


have 


Amplitude 


a(r, r) = a 0 


w n 


W[T — T„ 


-g w 2 (t — t w ) g 

I ^ 


r (t-t-uj) 


_ zl o o, ,u " —arctan T J w 

^ \ W 2 W 2 (T — T W ) Z 0 


(IV.21) 


Phase 


As can be seen, the amplitude falls off as e~ r2 /™ 2 and is clearly Gaussian while the 
phase term includes both r and r dependence. Equations (IV.20) and (IV.21) repre¬ 
sents the lowest-order transverse electromagnetic mode. Figure 16 shows the surface 


a i = 0 


a. i = 0 



100 


0 0 


0 0 


Figure 16. Real and imaginary parts of complex Gaussian mode at r = 0 with 
t w = 0.5. 

plots of the real and imaginary parts of Equation (IV.21). The real and imaginary 
parts of a(r, r) provides us with the amplitude and phase of the optical field since 

|a(r,r)| = ^{a} 2 + A{a} 2 (IV.22) 

and 

<j>(r, r) = arctan ^j-y (IV.23) 

where cj) is the optical phase. Burn patterns are recognizable through the power in the 
optical field, i.e., the optical field amplitude squared. In the following discussions we 
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will only plot the amplitude of the mode instead of the power (or amplitude squared) 
because we would lose detail in the plots. 


initial surface plot final surface plot 



0 0 x y 0 0 


Figure 17. Initial (r = 0) and final (r = 1) surface amplitudes for the lowest order 
mode where yellow represents highest largest amplitude. 



x 0 x lx 


Figure 18. Three panels shows cross sections of a Gaussian mode as it propagates 
from r = 0 —» 1. z 0 = 0.25 t w = 0. 


Figure 17 is a surface plot for a Gaussian mode at r = 0 and r = 1. Figure 
18 is composed of three panels. The first panel shows a cross section (x vs. y ) of the 
Gaussian mode at r = 0. The middle panel shows a cross section (y vs. r) as the 
mode propagates. The last panel shows the cross section of the mode at r = 1 (x vs. 
y). Notice how a small Rayleigh length has spread the mode considerably. 
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C. HIGHER ORDER MODES 


In order to account for more complicated variations along transverse directions, 
we can generalize the trial solution Equation (IV.3): 


a{r,r) = a Q g 


x 


W(T 


h 


WIT, 




zoq (t) 


(IV.24) 


The functions g and h will allow for a more intricate variance along transverse direc- 
tions. Again, we demand the same initial conditions, namely, a( r, 0) = a Q e . This 
implies, p(0) = 0, q( 0) = 1, and g(oc/w 0 ) = h(y/w 0 ) = 1 since the initial Gaussian 
form already expressed in the factor As before, our job now will be to find the 

unknown functions, p(r), q(r), g(x/w(r )), and h(y/w(r)) at later times. But before 
we begin taking derivatives, let us consider derivatives of g (and by analogy h) to 
illustrate the chain rule: 


and 


dg dg 8 / x \ 

dx 9(T) doc \w J 


1 dg 



d 2 g _ d dgi 
dx 2 dx dx 



dg dg d / x \ 

dr 9(4) dr \wJ 

_ ~ dg f 1 \ dib 


(IV. 25) 


(IV. 26) 


(IV. 27) 


As you can see, the primes over g denote a derivative with respect to its argument. 
This little fact will become handy later. Taking derivatives of Equation (IV.24) results 


in: 


38 



e 


(IV. 28) 


d 2 a 
dx 2 



Ag'hx 


wz 0 q 


2 gh Aghx 2 
z 0 q z 2 q 2 


-(p(r) + - 


■oQ(t) > 


d 2 a 

dy 2 



4 gh'y 


wz 0 q 


2 gh Aghy 2 
z 0 q z 2 q 2 


-(p(t)+ 


ioq(r) 


(IV. 29) 



1 ° 
— g'h w x 
w 2 


1 O o 

zz^gti w y-ghP -+ 
ur 


ghx 2 

z 0 q 2 


Q -f 


ghy 2 
z 0 q 2 


Q 


e 


(p(t)+ 


2 

- oq(T) 


(IV.30) 


Substituting Equations (IV.28), (IV.29), and (IV.30) into (IV.2), and dividing 
through by a 0 ghe we have the following unpleasant monster: 


1 g Ax g A , ill ~g 1 h 

zzy -rv- Ai—x— + -- - 

w 2 g wz 0 q g 

S_____^ 


w _ q' 1 h" Ay h' ,. w _ h' 

— r— -I --- Ai —v -1- 

w 2 g ih 2 h wz 0 q h w 2 h 


x only 


V only 


+ r 


z 2 q 2 


+ 


Ai 

z 0 Q 2 


- - AiP =0 (IV.31) 

z 0 q 

^-V-^ 


f 2 only 


no r dependence 


In order for Equation (IV.31) to have a solution at a given time r (i.e. a transverse 
plane), the x dependent terms should equal a constant, let’s say, —a. The y dependent 
terms should equal another constant, let’s say —15. The r 2 dependent terms should 
be identically zero since it must be valid for at any location in the transverse plane. 
Finally, the terms that don’t depend on r, should equal a + /5, so that the left hand 
side sums to zero. Solving the f equation, we have 


z 2 q 2 


+ 


Ai 


z a q* 

Af 2 r l ° 
— + iQ 


z 0 q Lz 0 


= 0 

= 0 


(IV.32) 
(IV.33) 


For this equation to hold for any f, the term in bracket must be identically zero. But 
we have seen this requirement before in the fundamental mode derivation. Applying 


39 



the initial conditions, we find that 


/ \ lT 

q( T ) = i + — 

Zq 


(IV.34) 


Returning to the x dependence of Equation (IV.31), we multiply through by 
gw 2 and collect terms in derivatives of g, giving us 


g" — 44: (-(- i w ) g + aw 2 g = 0 

V -o'l ) 


(IV.35) 


Let us evaluate the parenthetical term, recalling Equations (IV.17), (IV.18), and 

o 

noting that u: r /w 2 w, we find 


So now we have 


w . ° \ 

-b l W 

w 

~ 9 

(/- 

z oQ J 

W 2 0 

\w z 


1 



w 


// 

9 

4x , _ 2 

- —g + aw g 


(IV.36) 


(IV.37) 


The above equation is close to the form of the Hermite equation [12], y" — 2xy' + Xy = 
0. Recall that the primes denote a derivative with respect to the argument of the 
function. So a judicious change in variables will allow us to recast Equation (IV.37) 
in precisely the form of Hermite’s equation. Let 


a/24 _ 

= —i- > x = 

w a/2 


a? = \/2£>(-) 

w 


Using Equations (IV.38) and (IV.39) on (IV.37), results in 

dg 2 . dg aw 2 

w- 9 M + ~ 3=0 


(IV.38) 
(IV.39) 

(IV.40) 


A power series solution will solve Equation (IV.40) uniquely with initial conditions. 
But for the series to terminate, we must require aw 2 /2 to be equal to a nonnegative 
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even integer. With this requirement, the power series solution terminates into the 
well known Hermite polynomials: 


9«) = iUf) 





where m 


0 , 1 , 2 , 3 ... 


(IV.41) 


There will be a similar solution for the y terms of Equation (IV.31), i.e., h(r)) = 
H n (rj ) = H n (s/2y/w) where n = 0,1, 2, 3... The first ten Hermite polynomials are 
listed in Table II. 


H 0 (x) = 1 

Hi(x) = 2x 
H 2 (x) = 4x 2 - 2 

H 3 (x) = 8x 3 — 12a; 

H±(x) = 16a; 4 — 48 .t 2 + 12 

H 5 (x) = 32a; 5 — 160a’ 3 + 120a; 

H 6 (x) = 64a; 6 - 480a’ 4 + 720a; 2 - 120 

H r (x) = 128a’ 7 - 1344a; 5 + 3360a; 3 - 1680a; 

H 8 (x) = 256a; 8 - 3584a’ 6 + 13440a’ 4 - 13440a; 2 + 1680 

H 9 (x) = 512a’ 9 - 9216a’ 7 + 48384a’ 5 - 80640a; 3 + 30240a; 


Table II. Hermite Polynomials from [12] 


Returning to Equation (IV.31), we are left with finding a solution for p(r) 
For Hermite solutions, a and ft are constrained 


ailr 4 m 

= 2m —> a = —- 


2 

ftw 2 


ur 


„ . 4n 

= 2 n —> ft = 


2 w' 

o 

Solving for P. with the constraints, and Equations (IV.34) and (IV.17), 


-— - 4 ip= a + ft 
z 0 q 


(IV.42) 


P = (m + n) + -J— 
w z a q 

= i(m + n + 1) ° 


z 2 0 + T 1 


+ 


T 


5 2 + T 1 


(IV.43) 
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Integrating and applying initial conditions for p 


p(r) 



T 

i(m + n + 1) arctan — + In 
Z 0 



(IV.44) 


So we now have all the pieces to Equation (IV.24). Putting it altogether, Equations 
(IV.34), (IV.41), and (IV.44) into (IV.24), and after a little tidying up, we arrive 
finally at an expression for the (m,n) mode for the laser optical field: 


Amplitude 

_ 


a m Jr,T) = a 0 ^H m ( ^ ) H n ( ^ ) e~& 

W \w \w --- 

' ' ' Phase 


(IV.45) 


where the phase, 4> m , n , is given by 


(r,r) = 


f 2 r 


T 


— (m + n + 1) arctan — 


(IV.46) 


Figure 19 shows the “burn patterns” for various modes. 



Figure 19. |aoo|, |aoi|» |« 30 1» l a 3 i| (Left to right; Top to bottom) 
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D. NORMALIZATION 


If we attempt to superpose modes with high values of m or n, the Hermite 
polynomials with Gaussian envelopes for these higher modes have large coefficients 
as shown in Table II. To place all the modes on an equal footing we will normalize 
each of them to the power in the fundamental which is given by 


Po,o — 


oo poo 


ciofi\ 2 dxdy 


OO «/ —OO 


Substituting Equation (IV.20) we have 


P ~n 2W ° 
P °° - “°€^ 


2 r 0° 


e 


dx 


OO .2 

2» 

e ™ 2 1 


(IV.47) 


(IV.48) 


so then 

A, o = \a 2 c wl (IV.49) 

The power in an arbitrary mode is 


P = 

1 m.n 


oo poo 


Um,n dxdy 


(IV.50) 


oo •) — oo 


Substituting Equation (IV.45) we have 


P = 

1 m,n 


''OO 2 _ 2 

al^H^H^^dxdy 

-oo 


(IV.51) 


To evaluate this integral, we make use of the fact that Hermite polynomials are 
orthogonal on the range(—oo,oo) with respect to a weighting function [12], 





(IV.52) 


where the <) UJ is the Kronecker delta. To properly make use of the orthogonality 
condition, Equation (IV.52), we must first massage Equation (IV.51) into the correct 
form. Consider the factor 



(IV. 53) 
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(IV. 54) 


If we make the following change of variables, let 


\/2x _ aw 

a = - —} x = —-pz 

w V2 


so that, 


dx = —=da 

V2 


(IV. 55) 


Then Equation (IV.53) becomes 


w 


I H m (a) 2 e~ a da 
V 2 J — oo 


(IV.56) 

Employing the orthogonality condition, Equation (IV.52) when m! = m, we have 

(IV. 57) 

and similarly for the y factors. The power in the (m,n) mode then becomes 


w 


—V J H m (a) 2 e a da = il \ I ^2 in m\ 


P m ,n = | a 2 0 w 2 0 2 m 2 n m\n\ 
= P 0 , 0 2 m 2 n mM 


(IV. 58) 
(IV.59) 


So the power in the (m,n) mode can be normalized to the power in the fundamental 
by dividing it through by 2 m 2 n m\n\. A normalized expression for the optical field in 
the (m,n) mode is then 

Am,,n / , , (IV.60) 

V2 m 2 n m\n\ 

This expression places all the higher-order modes on an equal footing; all modes have 
the same power as the fundamental. 


E. FINDING COEFFICIENTS 

The total optical field is a superposition of all the separate normalized com¬ 
ponents, i.e., 

OO 

At{ r, t) = ^2 C m,nA m , n (r, r) (IV.61) 

m,n 

where C' m . n is the normalized coefficient for the (m,n) mode. Each term, or com¬ 
ponent, is a solution to the wave equation. So, the sum of terms is also a solution. 
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Moreover, the Hermite-Gaussian modes form a complete orthogonal basis set [16] so 
that any laser field satisfying the wave equation can be accurately represented by 
Equation (IV.61) with appropriate coefficients C m ^ n . Our task now is to determine 
the unknown normalized coefficients C' m . n . As with a Fourier series decomposition 

of a function, we can employ the Hermite analog of Fourier’s Trick 2 . To determine 

* 

the coefficients for an optical field, At , we multiply Equation (IV.61) by A m ^ n t and 
integrate over all x and y giving us 

/ OO /■•OO .. /‘OO /*oo °° .. 

I At dxdy / / ^ ] C m ^ n A m n A m '^ n i dxdy (I\ .62) 

oo J—oo J — ooJ—oo mn 

where the asterisk denotes taking the complex conjugate. Expanding out the right 
hand side, we have 


oo 1 

E r 1 w ° u u ~ 

V2 m 2 n m\n\ w 


1 ll, ° II IT - 

■ _ rim! 6 

V2 m 2 n mlnl w 


(IV.63) 

Notice the term that involves the product of the phases from the unprimed and primed 
optical field. Recalling Equation (IV.46), we then have 


•i ! . n ! — pi ( 0m, n -I ni 1 . n’ ) 


_ m.) + ;//.' //ii atrtan J j 


(IV.64) 


This equation does not depend on x or y, so we can pull it out of the integral. Re¬ 
turning to Equation (IV.63) and employing the orthogonality of Hermite polynomials 
when m' = m and n' = n, we have 


r»00 pOO 


At A m ji dxdy C n 


-oo J — oo 


2 TO 2 n m\n\ w 2 


n 2 'll] 2 

0 0 'wJ-2 m rn\ 


fiA/Vn! 


(IV.65) 


Combining this result with the left-hand side of Equation (IV.62), we have 


r 


oo poo 


7 T(l 0 U ' 0 •) — oo J — oo 


A t A m „ dxdy 


(IV.66) 


2 [14], Chapter 3 
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This equation is an expression for the normalized coefficients for our optical field A T . 
Our next task will be to write computer code to determine the coefficients of any 
waveform. 
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V. 


SIMULATIONS 


In this chapter we will describe how numerical simulation can solve the non¬ 
linear differential equations that have been developed. Non-linear differential equa¬ 
tions are exceedingly difficult to solve analytically but are relatively “easy” with a 
computer. 


A. MODELING THE FREE ELECTRON LASER 

In order to model the FEL we must simultaneously solve two coupled differ¬ 
ential equations: the pendulum equation, Equation (III.32), 

oo 

C= M cos(c + <j>) 


and the FEL wave equation, Equation (III.64), 


i - 


d 

+ d’T 


a(r,r) =< -je 


-K 


> 


The pendulum equation updates electron phases and phase velocities and the wave 
equation propagates the light self-consistently with the electron beam interaction. 
The Naval Postgraduate School’s Directed Energy and Electric Weapons Center uses 
an Apple Xserve cluster to model the FEL (see Figure 20). The cluster is specifically 
designed to model the FEL in four dimensions as well as model FEL performance due 
to perturbations or test validity of new designs. 

In order to numerically solve Equation (III.64), fast Fourier transforms are 
used in the code [17]. An output for a recent simulation is given in Figure 21. 
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Figure 20. NPS Apple Xserve cluster with 64 nodes / 128 processors. 
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nx=800, nt=BO, np=2B000, tmir=19, Wp=20 Fri May 26 IB:33:46 2006 

FEL WAVEFRONT EVOLUTICN 


|a(x,y,10)| 



j=60.0 a =0.22 
v 0 =12.B c£=0.22 
oo„=l. 2 0 ^= 0 .31 

i|!=O.S g*=0.31 

/=0 c =0.26 

;° =o <j- ve =1.4 

yo vy 

4=10 i^O.B 

z =0.42 r. =0.0 
o n 

r =40 r =9.S 
c 

e=3e-07 

0^=9 N=30 

A V r c =-O.BOO%( C ) 


P(n) 1.9xlO J 

G(n) (U 0.31 


0 ^ 


147 

2| 

4>(0,y,-9) 

|a 

* (0,y,10) 

Y 



. J 

] 

147 


0) 

(k) 

© 


-147 x 147=3 3 -147 x 147 -12 12 


P=1.66e+03 (2%), G=0.14B (97\), ^=0.0144 (4%), av=28.B 


Gmn=0.0279, Gmx=0.311, Gf=0.14B 
F=0.31 # G(th)=2, G(ext)=O.OS8, Af/'(=0.07B7 
v0=0.81 (0.6B), wl=14 (15), v2=14 (IB) 

M 2 =0.9B, c 2 (2,0)=0.31» c 2 (0,0)=0.31, c 2 (4,0)=0.27 
wbins=67, ebins=32 



Figure 21. FEL simulation output, (a) Optical mode evolution over 500 passes at 
right mirror (b) Surface plot of the optical field at the right mirror after last pass; 
(c) User defined parameters (see Table (III)); (d) Cross section of optical mode and 
electron beam as it propagates from r = 0 —> 1 on last pass; (e) Evolution of electron 
beam distribution in phase velocity over 500 passes; (f) Electron phase space at r = 1 
after last pass; (g) Optical power evolution over 500 passes; (h) Optical gain evolution 
over 500 passes; (i) Burn pattern of optical mode at left mirror after last pass; (j) 
Optical phase for at left mirror after last pass; (k) Burn pattern of optical mode at 
right mirror after last pass; (1) Optical phase for at right mirror for last pass; (m) 
Hermite-Gaussian coefficients for right mirror (Code is developed in this chapter) 
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Symbol 

Name 

Equation 

j 

Electron current density 

j = S7r 2 Ne 2 K 2 L 2 p/ r f 3 rnc 2 

v 0 

Electron phase velocity 

v 0 = [(k + k 0 )/3 z - k]L 

up 

Betatron frequency 

co is = 2 ttNK/'j 

T P 

Electron beam waist position 

t/ 3 = Zp/L 

Vo 

Electron beam shift 

y 0 = Y 0 y]n/L\ 

@yo 

Electron beam tilt 

0y O = Q\/L7r/\ 

a 

Initial optical field amplitude 

\a\ = AnNeK\E\L/^ 2 mc 2 

Z 0 

Rayleigh length 

z 0 Z 0 / L 

1 m 

Dimensionless mirror radius 

\J TC j\L j 

e 

Mirror edge loss 

See [18] 

Qn 

Cavity quality factor 

Q n = 1/(loss per pass) 


Electron beam radius 

( 7 X A e \J 7T jL\ 

&9y 

Electron beam angular spread 


<Jvc 

Spread in electron phase velocity 



due to emittance 

See [18] 

7 

Spread in electron phase velocity 



due to energy spread 

See [18] 


Optical mode waist position 

t w Z w /L 

Th 

Right mirror hole radius 

Th = Rh\/n/\L 

r c 

Mirror radius of curvature 

> 

II 

o 

5- 


Table III. User-defined FEL simulation parameters 


Each node of the Apple cluster can run a simulation such as the one shown 
above at a different values of u 0 simultaneously. We can then determine gain and 
extraction values at the optimum value of v 0 . This simulates what the FEL would do 
by mode competition. 

The simulation also provides us with an array describing the optical field, 
a(x,y). As mentioned in the previous chapter, the optical field is a superposition of 
Hermite-Gaussian basis set with appropriately chosen coefficients. Mirror distortion 
or electron beam shifts and tilts can cause the optical field to develop higher-order 
modes. We will now analyze these modes by finding their coefficients numerically. 
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B. FINDING HERMITE-GAUSSIAN COEFFICIENTS NU¬ 
MERICALLY 

Equation (IV.66) is the analytical expression for the normalized Hermite- 

Gaussian (HG) coefficients. We need to recast this equation in a form that is 

* 

computer-friendly. Since At and A m . n are complex, we can write: 


A t = . 1) + i.Y r 

(V.l) 

4* — 4 ! ' - } 4* 

s*m.n L ^m,n 

(V.2) 


where A^ and A' r are the real and imaginary parts of the total normalized optical 
field and A r m n and A' m n are the real and imaginary parts of the normalized (m,n) 
mode. Equation (IV.66) becomes 

2 


n — 

' m rt 





(At + iA! r)( A m,n ~ iA l mn )dx 


(V.3) 


Expanding out the factors and separating the coefficients into their real and imaginary 
parts, we have 


and 


CL 


CL 


2 

na 2 0 w 2 0 

2 

7T a 2 w 2 



+ a t A'm.n'l di dy 

- A 't A mj di d v 


(V.4) 

(V.5) 


|C„.„| 2 = (C;„,J 2 + (C” m ,J 2 (V.6) 

So, let us first describe the procedure for developing a numerical modal analysis tool 
in the C programming language: 


1. Construct Hermite-Gaussian basis arrays. 

2. Read in or form a user-defined optical field. 

3. Determine HG coefficients via Equations (V.4), (V.5), and (V.6) 
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Our expression for the optical field, Equation (IV.61) is a sum of the normalized basis 
at some time r. Recall that the optical field is a complex number so we must keep 
track of the real and imaginary parts. These analytical expressions are continuous in 
x and y so we recast the expression for the (m,n) mode into their real and imaginary 
parts over the discrete sites in the array. Therefore, in pseudo-code we have 


A r [m\[n\[x j][y k \ = 


. '^rr 
~ 1 1 771 

W 


n ^LfT 
Mo ~ 11 m, 

W 


y/2x. 


Hr, 


\/2 y k \ -% 

e * 2 cos 


w 


w 


r _M[_ 

WgW 2 


— (m + n + 1) arctan — ) (V.7) 


A l [m][n\[xj][y k ] = 


■y/2Xn 


Hr, 


w 


' V% Vk 
w 


JIk 

e * 2 


sin 


T XkT 

lUgW 2 


— (m + n + 1) arctan — ) (V.8) 


Z 0 


where A 1 and A r are four-dimensional arrays indexed by the mode number and posi¬ 
tion, given by (m,n) and (x 3 , y k ), respectively. We calculate the basis set for the first 
several modes and store these values in their arrays. We then read in the optical field 
that is passed as an array from the main FEL simulation or form a user-defined field 1 
by a linear combination of normalized Hermite-Gaussian modes using the arrays that 
been calculated. 

We determine the HG coefficients via Equations (V.4), (V.5), and (V.6). We 
approximate the double integral by nested loops over m, n, and over the number of 
array sites, nx. In pseudo-code, 

^ nx nx 

c r [m][n\ = —Y ( A T^j}[yk] + A T[xj\[yk] 

7T ciiwi ' *—r 
° ° k j 

(V.9) 

and 


nx nx 


C*[m][n] = — bEEW pj\\yk\ A r [m][n\[x j ][y k ] - A T [xj][y k \ ^‘[m][n][ij][y fc ]) 

(V.10) 


7ia 2 0 w 2 0 k 


b4 user-defined field can be used as a check on the code to ensure that it is operating properly. 
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The square of the magnitude of the HG coefficient for the (m,n) mode is then given 
by Equation (V. 6 ), that is, 

C[m][n ] 2 = (C r [m] [n ]) 2 + (C' ,: [m][n ]) 2 (V.ll) 


It is convenient to normalize the square of the magnitudes to the power of the 
total field A t . The power in the total field is given by 


Pt = 


poo 

8__ 

/—00 J 

— OO 


(V-12) 


Substituting Equation (IV.61) and working out the integrals we find that 

— 00 

Pt = \4alY. C ™.n ( v '!3) 

m,n 

Recall that the power in the fundamental mode is given by Po,o = 7 ra 2 t (; 2 /2 so that 
the total power in the optical field is then (in terms of the fundamental) 


Pt = -Po,o Cl., 

m,n 

The fractional power spectrum, S[m][n], is then 


S [m] [n = 


C[m] [n ] 2 

WM 


(V.14) 


(V.15) 


where ■ C\i] [j] 2 is the sum of all the squares of the normalized mode coefficients. 

To illustrate how the code works, consider the following examples illustrated in 
Figures 22 and 23. In Figure 22 (a), the diagnostic correctly identifies the fractional 
power spectrum as 50% in both y4 0i2 and % 2 ,0 modes. Incidentally, the superposition 
of equally weighted % 0 ,2 and % 2i0 , with no relative phase difference between them, 
results in a pure L\ Laguerre-Gaussian mode 2 . 

In Figure 22 (b), the diagnostic correctly identifies the fractional power spec¬ 
trum for the superposition of % 0 , 0 i ^4i,i, and h % 3 with C^o = 1, Cq 1 = 2, and C 3i 3 = 3 


2 To see this: ask Professor Colson and he will buy you a beer. 
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as 7%, 28%, and 64%, respectively. A screen capture of Cmn.c is shown in Figure 
23 for this case. As you can see, the program outputs the Rayleigh length and mode 
waist. It then displays the mode radius at r = 0, the power in the fundamental, and 
total power in the arbitrary field. The fractional power spectrum is then calculated 
and displayed. The code then propagates the optical field in free-space and performs 
the same analysis on the final waveform. 

In Figure 22 (c), the diagnostic correctly identifies the fractional power spec¬ 
trum for the superposition of equally weighted A 0 ,0) A 0i2 , and A 0 ,4 modes. We will see 
a similar combinations of even ‘nT and £ n’ valued modes when we investigate mirror 
distortions. 

The diagnostic developed for this thesis for Hermite-Gaussian modal analysis 
was then incorporated into the main NPS FEL simulation code as a function call. 
In the next chapter, we will see this diagnostic in action when we simulate mirror 
distortions. 
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(a) 



• II 




Aq 2 + Ao 


Code results 

50% Aq 2 
50% A 20 



7% A 00 
29% A n 
64% A 33 


33% Aq 0 
33% Aq 2 
33% Aq 4 


Figure 22. Graphics from Cmn.c output was drawn with MatLab (a) Superposition 
of equally weighted 4 0 ,2 and A 2 ,o- (b) Superposition of A 0j o, 4.^1, and 4 3i3 with 
C 0 ,o = 1, CV, = 2, and C 3i3 = 3. (c) Superposition of equally weighted 4 0 0 , 4 0i2 , 
and 4 0 4 . 
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Files opened 

Initializalions 

zo = 0.5 tanw = 0. 5 

Loading np basis sets at tan = 0. 

w(0) = 1 

P 00(tan = 0) = 78.5398 
P_T(0) = 1099.56 

Outputting initial wavefront to file 
Calculating initial Crms 

Normalized initial cmn amplitudes squared 
Cmn2[0][0] = 0.0711286 
Cmn2[1][1] = 0.285714 
Cmn2[3][3] = 0.642857 

Propagating waveform. . . 

w(l) = 1 

P_00(1) = 78.5398 
P_T(1) = 1119.46 

Nomalized initial cmn amplitudes squared 
Cmn2[0][0] = 0.0702047 
Cmn2[1][1] = 0.28182 
Cmn2[3][3] = 0.639894 


Closing files... 

spanl05-250:~/Desktop/Cnrrent Code rvigilpbg4$ 


Figure 23. Output from diagnostic tool Cmn.c 
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VI. 


MIRROR DISTORTIONS 


A. COLD CAVITY THEORY 

In the pursuit of a weapons class FEL, high power and high intensity can lead 
to mirror distortions if heat deposited by the laser light to the mirror substrate isn’t 
removed by cooling in a uniform manner. Distortions then cause higher-order modes 
to appear. It is important to quantify the fractional power spectrum among the 
higher-order modes because doing so uniquely identifies the laser field and may allow 
for corrective optics to single out the best modes for FEL lethality. We will investigate 
three different cases of mirror distortions on the outcoupling mirror: hyperbolic, 
ellipsoidal, and spherical. Hyperbolic and ellipsoidal distortions can be classified as 
a mirror astigmatism. It occurs when the radius of curvature of the spherical mirror 
along an axis differs by some small amount when compared to the perpendicular axis. 
For simplicity we chose these axes to line up with x and y. Spherical distortion is 
not astigmatism because both axis are changed in exactly the same amount and in 
the same direction. Spherical distortions do change the Rayleigh length. Before we 
consider the effects of mirror distortions on FEL operation, let us explore cold-cavity 
effects (i.e., no electron beam present) of mirror “errors” to get a flavor for how the 
FEL might respond. Keep in mind that cold-cavitv theory discussed in this section 
keeps both mirrors the same. The dimensionless radius of curvature of the two ideal 
mirrors determines the Rayleigh length, z 0 [18], and is given bv 1 


~ = T m 2z 0 
C 2 + Tm 


(VL1) 


where r m is the dimensionless mirror separation. Factoring out r m /2 and defining a 
small parameter k. = we then have 


Rra / -i , \ 

= T d + 0 


(VI.2) 


Symmetric confocal cavity. See [19], chapter 12. 
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where k « 1 for r m >> z 0 . Typically, r m fs 10 and z 0 fa 0.1, so k fa 4 x 1CT 4 . What 
if an actual mirror differs from the ideal mirror by a small amount? This could occur 
because mirror heating has distorted both mirrors, for instance. We can quantify this 
distortion as follows 

r c -> f c (l + e) (VI.3) 

where e is a small deviation from the original radius of curvature r c . Combining these 
relations we have 

f c = ^(1 + «)(1 + C ) (VI-4) 

r c = ^(1 + k + e + kc) (VI.5) 

Since k and c are small, we can approximate this equation to first order giving us 

r c -^-(1 + k, + e) (VI.6) 

We now ask, “What is going on with the mode radius at the waist and at the 
outcoupling mirror with respect to k and r TO ?” The waist is given by 



We can solve Equation (VI. 1) for z 0 giving us 



and substitute that into our expression for the waist giving us 



Substituting Equation (VI.6), we now have 


K (t) , (|;y ( 1+b+£) “ 1 ) 

W 0 « ^ — J (K, + e) 4 


(VI.7) 

(VI.8) 

(VI.9) 

(VI.10) 

(VI.11) 
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Let us now consider the mode radius at the mirror. The mode radius is given 
by Equation (IV.17). At the outcoupling mirror, r = r m / 2, we then have 


w(r m / 2) = w. 



(VI.12) 


Substituting Equations (VI.6) and (VI.8) in for z 0 and after a little algebra we arrive 
at 


Wm = tPf (VI.13) 

V« + £ 

In both Equations (VI.11) and (VI.13), there is a problem when the mirror distortion 
e —> —k. Let us plot these two equations with respect to e/n. From Figure 24 we 


Waist radius Mode radius at mirror 



(a) (b) 


Figure 24. Cold-cavity, (a) Waist radius variation (b) Mode radius variation at the 
outcoupling mirror. 

can see that if the error, e is less than —k the waist radius and mode radius becomes 
imaginary. If e = —k, the waist radius becomes zero and the mode radius at the 
mirror becomes infinity. So, as far as cold-cavity analysis predicts, mode radius can 
be sensitive to mirror distortions that are on the order of —k. So we may expect 
actual FEL operation to follow a similar trend. 
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We chose to simulate the effects of mirror distortions on two different .JLab 
oscillator configurations: Electromagnetic (EM) wiggler and STI Optronics wiggler 
designs. The following table highlights their major differences. 


EM wiggler STI wiggler 


Undulator period, A 0 (cm) 

8 

5.5 

Undulator parameter, K 

0.5 

1.36 

Undulator length, L (m) 

2.3 

1.65 

No. Periods, N 

29 

30 

Mirror Separation, S (m) 

30 

50 

Quality factor, Q n 

25 

9 

Ravleigh length. Z 0 (cm) 

90 

70 

K (%) 

0.35 

0.2 


Table IV. JLab EM and STI wiggler comparison 


B. HYPERBOLIC DISTORTION 

In the case of hyperbolic distortion, we vary the mirror distortion on the 
outcoupling mirror in the following way: 

A fey = - A r ci (VI. 14) 

where A f C y and Af c ^ are the changes of outcoupling mirror radius-of-curvature in the 
x and y directions, respectively. First let us consider the EM wiggler design. The 
results for gain and extraction are shown in Figure 25. Notice how both the gain and 
extraction graphs are symmetric about zero e as should be expected for a hyperbolic 
distortion. Although, cold cavity analysis predicts catastraphic results when mirror 
distortions are less than —k, we see that with FEL operation both the gain and 
extraction curves continues well beyond this cold cavity limit (—k is indicated by the 
red tick mark). 

Notice also a peak in gain and extraction as the magnitude of e is increased. 
It seems counterintuitive that increasing hyperbolic mirror distortion would improve 
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e (%) (= A r c /r c ) e (%) (= A r c /r c ) 

(a) (b) 

Figure 25. (a) Gain trend for hyperbolic distortion for the JLab EM wiggler. (b) 
Extraction trend for hyperbolic distortion for the JLab EM wiggler 

FEL operation. This can be explained if we look at at Figure 25 (b). If we compare 
the burn patterns for e = 0 and e = 0.5%, we can see that the mode shape at the 
mirror (yellow) is wider for the e = 0.5% case than for the e = 0% case. If the 
mode is more spread out at the mirror, we can conclude that the mode was more 
tightly focused in the undulator leading to a greater interaction between the optical 
mode and the electron beam leading to greater gain and extraction. This beneficial 
effect does not go on indefinitely. Eventually greater distortion will cause the mode 
to become so large that the mode will be clipped by the outcoupling mirror reducing 
gain and extraction. An example of a mode being clipped is shown in Figure 25 (a) 
for the e = — 2% case. We also show the modal coefficients that comprise this data 
point. Notice that the ‘n’ modes are weighted more heavily than the ‘nT modes since 
the burn pattern is spread out in the y direction and that only even modes are present 
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since the original wave shape is assumed to be a fundamental Gaussian. Only the 
even modes overlap the origin (x = 0, y = 0) where the electron beam is amplifying 
the light. 



Figure 26. (a) Gain trend for hyperbolic distortion for the .JLab STI wiggler. (b) 
Extraction trend for hyperbolic distortion for the JLab STI wiggler 

We next considered the STI wiggler design. The results for gain and extraction 
are shown in Figure 26. As with the EM wiggler case, we see that the gain and the 
extraction are symmetric. The gain is much larger in the STI wiggler design due to a 
larger value of the undulator parameter, K . The mirror separation is greater, thereby 
yielding a somewhat smaller value for k (. k E m = 0.35% vs. k ST i = 0.2%). We also see 
a rise in gain and extraction for distortions on the order of \k\ as in the EM wiggler 
case. Notice in Figure 26 (b) for e = 0.5%, how the coefficients are strongly weighted 
in the ‘m’ modes due to a spreading out in the x direction. Again, only the even 
modes are driven by the electron beam at the origin. 
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C. ELLIPSOIDAL DISTORTION 


In ellipsoidal distortion, we vary the mirror error on the outcoupling mirror in 
the following way: 

vary A r c§ , A f ci . = 0 (VI. 15) 

As before, let us consider the EM wiggler results. Notice in Figure 27, both the gain 



5 ^ 


Figure 27. (a) Gain trend for ellipsoidal distortion for the JLab EM wiggler. (b) 
Extraction trend for ellipsoidal distortion for the JLab EM wiggler 




and extraction curves are not symmetric, as should be expected for an asymmetric 
distortion. Both gain and extraction curves have the same general shape as Figure 
24 (a). Unlike the cold cavity case, the FEL operates well beyond e < —k. As 
e is decreased further, higher-order modes appear as shown in the coefficient chart. 
Conversely, as e is increased (e > —k), the mode remains near the fundamental mode. 

In Figure 28, we have the results for ellipsoidal distortions in the STI wiggler. 
Both curves are once again asymmetric. Notice how the gain curve reaches its max- 
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Figure 28. (a) Gain trend for ellipsoidal distortion for the JLab STI wiggler. (b) 
Extraction trend for ellipsoidal distortion for the JLab STI wiggler 

imum at e = —1.5% which is several times greater than cold-cavity theory predicts 
the resonator should fail. The extraction curve is not as smooth as the EM wiggler 
extraction curve. But the general shape is still present, i.e., maximum nigh —e and a 
relatively flat extraction in the positive e direction. 

D. SPHERICAL DISTORTION 

In spherical distortion, we vary the mirror distortion on the outcoupling mirror 
in the following way: 

vary A r C y = A r ci (VI. 16) 

Spherical distortion is not a misnomer even though the mirror remains cylindrically 
symmetric, because the highly reflecting mirror is not altered. The high power laser 
beam is assumed to have caused a symmetric distortion, as in Equation (VI.16). Res- 


64 



















onator cavity instabilities are introduced as the distortion is increased or decreased. 
Notice in Figure 29, the burn patterns are cylindrically symmetric. Again, we see 



Figure 29. (a) Gain trend for spherical distortion for the .JLab STI wiggler. (b) 
Extraction trend for spherical distortion for the .JLab STI wiggler 


good gain and extraction well into the regime where cold-cavity theory predicts a 
non-functioning FEL. Modal decomposition reveals an equal weighting of lower-order 
modes yielding burn patterns with central holes in them. Spherical distortion varies 
the effective cavity Rayleigh length. In Figure 30, we notice the burn patterns caused 
by varying Rayleigh length are cylindrically symmetric and comprised of equally 
weighted modes similar to spherical distortions. In Figure 30 (a), gain increases as 
we go to shorter Rayleigh lengths because the mode is more focused in the interac¬ 
tion region (see Figure 14) increasing the value of the filling factor. Since, gain is 
proportional to the filling factor, gain goes up. Shorter Rayleigh lengths expands the 
mode on the mirrors. In Figure 30 (b), the green curve shows a decreasing intensity 
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Figure 30. (a) Gain trend for Rayleigh length variation in the JLab STI wiggler. (b) 
Extraction and mirror intensity trends for Rayleigh length variation in the JLab STI 
wiggler 


with decreasing Rayleigh length. This is good for FEL optics because the decreased 
intensity will mitigate mirror distortions caused by uneven heating or cooling. 
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VII. 


CONCLUSION 


In this thesis, we have reviewed free electron laser theory and for the first 
time developed analytical solutions to quantify Hermite-Gaussian higher-order modes, 
developed a diagnostic for modal analysis, and investigated tolerance limits on mirror 
distortions. Hermite-Gaussian modes form a complete and orthogonal basis set that 
can fully describe any optical field and is a natural basis set when describing mirror 
distortions that do not possess cylindrical symmetry. 

Modal analysis is a useful tool for describing mode quality. We are able to 
identify the fractional power spectrum in higher-order modes. This information may 
be used by an adaptive optics system to single out the best modes for propagation 
and lethality. 

Tolerance limits on mirror distortions are several times greater than what cold 
cavity theory predicts. This is due to the electron beam not allowing the optical mode 
to shrink to zero size. The electron beam stabilizes the operation of the FEL. Sim¬ 
ulations indicate JLab FEL should tolerate mirror distortions several percent, which 
is much greater than cold-cavity predictions. Finally, we recommend the following 
future work: 

• Investigate Laguerre-Gaussian modes as another complete and orthogonal ba¬ 
sis set to include azimuthal dependence. 

• Incorporate both Hermite-Gaussian and Laguerre-Gaussian modal decompo¬ 
sition in the complete four dimensional FEL simulation program. 

• Utilizing the diagnostic tool, investigate higher-order modes that appear due 
to mirror shifts/tilts and electron beam shifts/tilts. 

• Choose a catchy and snappy name for NPS FEL simulation software....something 
other than wavens. 
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APPENDIX A. CODES 


We present the following programs used for finding Hermite-Gaussian coef¬ 
ficients of an arbitrary optical field and the MatLab scripts used in plotting. Line 
numbers are included for convenience. 


1. CMN.C 

/ * start of Cmn.c */ 

1 // LCDR Ricardo Vigil 

2 // Thesis work 

3 // June 2006 graduation date 

4 // 13 major revisions 

5 

6 /* Program functions: 

7 - Combines up to ten user defined Hermite Gaussian modes 

8 to form an initial optical field 

9 - Decomposes the initial optical field into its constituent HG modes 

10 - Propagates the waveform in free space 

11 - Decomposes the propagated waveform into its constituent HG modes. 

12 

13 input files: Cmn.in 

14 

15 output files: 

16 initial.out - (x,y) initial amplitude in optical field, i.e, at tau = 0 

17 final.out - (x,y) final amplitude in optical field, i.e., at tau = 1 

18 ytau.out - (y,t) amplitude cross section in (y, tau = 0 to 1) 

19 

20 Matlab graphics file: Coeffs.m 

21 Creates 5 graphics 

22 1. Amplitude burn pattern at tau = 0 

23 2. Amplitude surface plot at tau = 0 

24 3. Amplitude burn pattern at tau = 1 

25 4. Amplitude surface plot at tau = 1 

26 5. Propagation cross section in y vs tau 

27 

28 */ 

29 

30 #include <stdio.h> 

31 #include <stdlib.h> 
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32 #include <math.h> 

33 

34 

35 /* these arrays are defined outside mainO due to segmentation faults */ 

36 

37 double Abinitialr[10][10][200][200]; // real: (m,n) tau = 0 

38 double Abinitiali[10][10][200][200]; // img: (m,n) tau = 0 

39 double Abfinalr[10][10][200][200]; // real: (m,n) tau = 1 

40 double Abfinali[10][10][200][200]; // img: (m,n) tau = 1 

41 

42 int main(){ 

43 

44 double P0=0.0, P1=0.0, rsq=0.0, dx=0.0, dy=0.0, C=0.0, wi=0.0, wf=0.0; 

45 double Ao=0.0, zo=0.0, wo=0.0, tau=0.0, tauw=0.0, dtau=0.0, phi=0.0; 

46 double c2sumlnitial=0.0, c2sumFinal=0.0, templ=0.0, temp2=0.0; 

47 double Hx=0.0, Hy=0.0, ampinitial=0.0, ampfinal=0.0, ampy=0.0, K=0.0; 

48 double pi=3.1415926, u=0.0, v=0.0; 

49 

50 double cinitialr[10][10]; // real: inital coeff for the mth nth mode 

51 double cinitiali[10][10]; // img: intial coeff for the mth nth mode 

52 double c2maginitial[10][10];// initial magnitude squared for (m,n) mode 

53 double cfinalr[10][10]; // real: final coeff for the mth nth mode 

54 double cfinali[10][10]; // img: final coeff for the mth nth mode 

55 double c2magfinal[10][10]; // final magnitude squared for (m,n) mode 

56 double iS[10][10]; // initial power spectrum ratio for (m,n) mode 

57 double fS[10][10]; // final power spectrum ratio for (m,n) mode 

58 

59 double Powerinitial[10] [10]; //Power in (m,n) normalized. 

60 double Powerfinal[10][10]; //Power in (m,n) normalized. 

61 

62 double Atr[200] [200]; //real: total real optical field 

63 double Ati [200] [200]; //img: total img optical field 

64 double Atr_old[200][200]; //old values: real 

65 double Ati_old[200][200]; //old values: img 

66 

67 int i=0,j=0,k=0,1=0,nx=0,W=0,x=0,y=0; 

68 double 00=0.0,01=0.0,02=0.0,03=0.0,04=0.0; 

69 double c5=0.0,c6=0.0,c7=0.0,c8=0.0,c9=0.0; 

70 int m0=0,ml=0,m2=0,m3=0,m4=0,m5=0,m6=0,m7=0,m8=0,m9=0; 

71 int n0=0,nl=0,n2=0,n3=0,n4=0,n5=0,n6=0,n7=0,n8=0,n9=0,count1=0; 

72 unsigned int f; 

73 unsigned int factorial(unsigned int a); 
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74 

75 int m[10] , n[10] ; 

76 double c[10]; 

77 

78 FILE *input, ^initial, *final, *ytau; 

79 input = fopen("Cmn.in","r"); 

80 initial = fopenO'initial.out" , "w"); 

81 final = fopenO'f inal. out" , "w"); 

82 ytau = fopen("ytau.out","w"); 

83 

84 printf("Files opened\n"); 

85 

86 fscanf (input, "%lf %d °/ 0 d" ,&cO,&mO,&nO); 

87 fscanf(input,"%lf %d %d"j&clj&mlj&nl); 

88 fscanf(input,"%lf %d %d" J &c2,&m2,&n2); 

89 fscanf(input,"%lf %d %d",&c3,&m3,&n3); 

90 fscanf(input,"%lf %d %d",&c4,&m4,&n4); 

91 fscanf(input,"%lf %d %d",&c5,&m5,&n5); 

92 fscanf(input,"%lf °/ 0 d 0 / 0 d" ,&c6,&m6,&n6) ; 

93 fscanf (input, "%lf %d °/ 0 d" ,&c7,&m7,&n7); 

94 fscanf(input,"%lf °/ 0 d 0 / 0 d" ,&c8,&m8,&n8) ; 

95 fscanf (input, "%lf %d °/ 0 d" ,&c9,&m9,&n9); 

96 fscanf (input, "°/ 0 lf °/«d °/ 0 d" ,&Ao ,&¥,&nx); 

97 fscanf (input, "°/ 0 lf °/«lf %lf " ,&zo ,&tauw,&dtau); 

98 

99 /* initializations-*/ 

100 

101 printf("Initializations'^"); 

102 

103 c[0]=c0; c[1] =cl; c[2]=c2; c[3]=c3; c[4]=c4; 

104 c[5]=c5; c[6]=c6; c[7]=c7; c [8] =c8; c[9]=c9; 

105 

106 m[0]=m0; m[l]=ml; m[2]=m2; m[3]=m3; m[4]=m4; 

107 m[5]=m5; m[6]=m6; m[7]=m7; m[8]=m8; m[9]=m9; 

108 

109 n[0]=n0; n[l]=nl; n[2]=n2; n[3]=n3; n[4]=n4; 

110 n[5]=n5; n[6]=n6; n[7]=n7; n[8]=n8; n[9]=n9; 

111 

112 dx = (1.0*W)/(1.0*nx); 

113 C = dtau/(4.0*dx*dx); 

114 wo = sqrt(zo); 

115 printf("zo = %g\t tauw = %g\n", zo,tauw); 
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116 

117 for(k=0;k<10;k++) 

118 for(1=0;1<10;1++) { 

119 Powerinitial[k][1]=0.0; 

120 Powerfinal[k][1]=0.0; 

121 cinitialr[k][1] = 0.0; 

122 cinitiali[k][1] = 0.0; 

123 c2maginitial[k][1] = 0.0; 

124 cfinalr[k] [1] = 0.0; 

125 cfinali[k][1] =0.0; 

126 c2magfinal[k][1] = 0.0; 

127 iS[k] [1] = 0.0; 

128 fS[k] [1] =0.0; > 

129 

130 for(x=0;x<nx;x++) 

131 for(y=0;y<nx;y++) { 

132 Atr[x] [y] = 0.0; 

133 Ati [x] [y] = 0.0; 

134 Atr_old[x][y] = 0.0; 

135 Ati_old[x][y] = 0.0; } 

136 

137 for(k=0;k<10;k++) 

138 for(1=0;1<10;1++) 

139 for(i=0;i<nx;i++) 

140 for(j=0;j<nx;j++){ 

141 Abinitialr[k][1][i] [j] = 0.0; 

142 Abinitiali[k][1][i] [j] = 0.0; 

143 Abfinalr[k] [1] [i] [j] = 0.0; 

144 AbfinaliEk]El]Ei]Ej] = 0.0; } 

145 


146 /* Initializaitons complete -*/ 

147 /* Loading basis set-*/ 

148 

149 printf("\nLoading up basis sets at tau = 0.\n\n") 


150 

151 tau = 0.0; 

152 wi = sqrt( zo + (tau-tauw)*(tau-tauw)/zo ); 

153 printf ("w(7og) = %g\n" ,tau, wi) ; 

154 

155 for(k=0;k<10;k++) 

156 for(1=0;1<10;1++) 

157 for(x=0;x<nx;x++) 
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158 for(y=0;y<nx;y++){ 

159 u = (sqrt(2)/wi)*(x - nx/2 + 0.5) * dx; 

160 v = (sqrt(2)/wi)*(y - nx/2 + 0.5) * dx; 

161 if (k==0) Hx=l.0; 

162 if (k==l) Hx=2*u; 

163 if (k==2) Hx=4*pow(u,2)-2; 

164 if (k==3) Hx=8*pow(u,3)-12*u; 

165 if (k==4) Hx=16*pow(u,4)-48*pow(u,2)+12; 

166 if (k==5) Hx=32*pow(u,5)-160*pow(u,3)+120*u; 

167 if (k==6) Hx=64*pow(u,6)-480*pow(u,4)+720*pow(u,2)-120; 

168 if (k==7) Hx=128*pow(u,7)-1344*pow(u,5)+3360*pow(u,3)-1680*u; 

169 if (k==8) Hx=256*pow(u, 8)-3584*pow(u,6)+13440*pow(u,4)-13440* 

170 pow(u,2)+1680; 

171 if (k==9) Hx=512*pow(u,9)-9216*pow(u,7)+48384*pow(u,5)-80640* 

172 pow(u,3)+30240*u; 

173 if (1==0) Hy=l.0; 

174 if (1==1) Hy=2*v; 

175 if (1==2) Hy=4*pow(v,2)-2; 

176 if (1==3) Hy=8*pow(v,3)-12*v; 

177 if (1==4) Hy=16*pow(v,4)-48*pow(v,2)+12; 

178 if (1==5) Hy=32*pow(v,5)-160*pow(v,3)+120*v; 

179 if (1==6) Hy=64*pow(v,6)-480*pow(v,4)+720*pow(v,2)-120; 

180 if (1==7) Hy=128*pow(v,7)-1344*pow(v,5)+3360*pow(v,3)-1680*v; 

181 if (1==8) Hy=256*pow(v,8)-3584*pow(v,6)+13440*pow(v,4)-13440* 

182 pow(v,2)+1680; 

183 if (1==9) Hy=512*pow(v,9)-9216*pow(v,7)+48384*pow(v J 5)-80640* 

184 pow(v,3)+30240*v; 

185 

186 rsq = ((x-nx/2+0.5)*(x-nx/2+0.5)+(y-nx/2+0.5)*(y-nx/2+0.5))*dx*dx; 

187 phi = (rsq*(tau-tauw)/(wo*wo*wi*wi)-(k+1+1)*atan((tau-tauw)/zo)); 

188 

189 Abinitialr[k][1][x][y]=Ao*wo*Hx*Hy*exp(-rsq/(wi*wi))*cos(phi)/wi* 

190 sqrt(1/(pow(2,k)*pow(2,1)*factorial(k)*factorial(1))); 

191 

192 Abinitiali[k][1][x][y]=Ao*wo*Hx*Hy*exp(-rsq/(wi*wi))*sin(phi)/wi* 

193 sqrt(1/(pow(2,k)*pow(2,1)*factorial(k)*factorial(1))); 

194 

195 Powerinitial[k][1] = Powerinitial[k][1]+ 

196 Abinitialr[k][1][x][y]*Abinitialr[k] [1][x][y] + 

197 Abinitiali[k][1][x][y]*Abinitiali[k][1] [x] [y] ; } 

198 

199 /* This serves as a check to ensure modes add up in the expected way */ 
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200 

201 printf ("P_00(tau = 0) = °/ 0 g\n" ,Powerinitial [0] [0]); 

202 

203 

204 

205 /* forming initial At 0 tau = 0-*/ 

206 for(i=0;i<10;i++) 

207 for(x=0;x<nx;x++) 

208 for(y=0;y<nx;y++) { 

209 Atr[x] [y] = Atr[x] [y] + c [i] *Abinitialr [m[i] ] [n[i] ] [x] [y] ; 

210 Ati[x][y] = Ati[x][y] + c [i] *Abinitiali [m[i] ] [n[i] ] [x] [y] ; } 

211 

212 /* finding total power in At */ 

213 

214 for(x=0;x<nx;x++) 

215 for(y=0;y<nx;y++) 

216 P0 = P0 + Atr[x] [y]*Atr[x] [y] + Ati [x] [y] *Ati [x] [y] ; 

217 

218 printf ("P_T(0) = */.g\n",P0); 

219 printf("Outputting initial wavefront to file\n\n"); 

220 

221 /* Output At cross section to file and saving old values*/ 

222 

223 for(x=0;x<nx;x++) { 

224 for(y=0;y<nx;y++) { 

225 ampinitial = sqrt( Atr [x] [y] *Atr [x] [y] + Ati [x] [y] *Ati [x] [y] ); 

226 Atr_old[x] [y] = Atr[x] [y] ; 

227 Ati_old[x] [y] = Ati[x][y]; 

228 fprintf(initial, "%lf ",ampinitial); } 

229 fprintf(initial, "\n"); > 

230 

231 /* Cmn calculation for initial waveform-*/ 

232 

233 printf("Calculating initial Cmns\n"); 

234 

235 K = (2*dx*dx)/(pi*zo); 

236 for(k=0;k<10;k++) 

237 for(1=0;1<10;1++) 

238 for(x=0;x<nx;x++) 

239 for(y=0;y<nx;y++) { 

240 cinitialr[k] [1] = cinitialr[k] [1] + K * 

241 (Atr[x] [y] *Abinitialr[k] [1] [x] [y] + 
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242 Ati[x][y]*Abinitiali[k][1][x] [y]); 

243 

244 cinitiali[k][1] = cinitiali[k] [1] + K * 

245 (Ati[x][y]*Abinitialr[k][1][x] [y] - 

246 Atr[x][y]*Abinitiali[k][1][x][y]);> 

247 

248 

249 /* calculating Cmn squared */ 

250 for(k=0;k<10;k++) 

251 for(1=0;1<10;1++) { 

252 c2maginitial[k] [1] = (cinitialr[k][1]*cinitialr [k] [1] + 

253 cinitiali[k][1]*cinitiali[k][1]); 

254 c2sumlnitial = c2sumlnitial + c2maginitial[k] [1] ; > 

255 

256 printf("\nNormalized initial cmn amplitudes squared\n"); 

257 

258 for(k=0;k<10;k++) 

259 for(1=0;1<10;1++) { 

260 iS[k][l] = c2maginitial[k] [l]/c2sumlnitial; 

261 if (iS[k][l] > 0.01) printf ("Cmn2[%d] [°/ 0 d] = °/ 0 g\n" ,k,l,iS[k] [1]) ; } 

262 

263 

264 /* propagation of At from tau = 0 to tau= 1 */ 

265 

266 printf("\nPropagating waveform...\n\n"); 

267 

268 for(tau=0.0;tau<=0.9999;tau=tau+dtau) {/* start tau loop */ 

269 

270 /* wave equation over x,y and tau */ 

271 for(y=l;y<nx-l;y++) 

272 for(x=l;x<nx-l;x++){ 

273 Atr[x] [y] = Atr [x] [y] -C* (Ati_old[x+l] [y] +Ati_old[x-l] [y] + 

274 Ati_old[x][y+1]+Ati_old[x][y-l]-4.0*Ati_old[x][y]); 

275 

276 Ati[x][y] = Ati[x] [y]+C*(Atr_old[x+l] [y]+Atr_old[x-l] [y] + 

277 Atr_old[x][y+l]+Atr_old[x][y-l]-4.0*Atr_old[x] [y]); } 

278 

279 /* saving old values */ 

280 for(x=0;x<nx;x++) 

281 for(y=0;y<nx;y++){ 

282 Atr_old[x] [y] = Atr[x] [y] ; 

283 Ati_old[x] [y] = Ati[x][y]; } 
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284 

285 /* y vs. tau output */ 

286 if (countl==100) { 

287 count1=0; 

288 for(y=0;y<nx;y++) { 

289 ampy = sqrt (Atr[nx/2][y]*Atr[nx/2] [y] + Ati[nx/2][y]*Ati[nx/2][y]); 

290 fprintf (ytau, "°/ 0 g ",ampy); } } 

291 

292 fprintf(ytau,"\n"); 

293 

294 countl = countl +1; } /*end tau loop */ 

295 

296 /* output for final waveform */ 

297 for(x=0;x<nx;x++) { 

298 for(y=0;y<nx;y++){ 

299 ampfinal = sqrt( Atr[x] [y]*Atr[x] [y] + Ati [x] [y] *Ati [x] [y] ); 

300 fprintf(final,"%g ",ampfinal); } 

301 

302 fprintf(final,"\n"); > 

303 

304 

305 /* Calculating cmns for final waveform (i.e., after propagation)- */ 

306 /* Must create basis set at tau = 1.0 */ 

307 

308 wf = sqrt( zo + (tau-tauw)*(tau-tauw)/zo ); 

309 printf ("w(°/ 0 g) = °/ 0 g\n" ,tau, wf) ; 

310 

311 

312 for(k=0;k<10;k++) 

313 for(1=0;1<10;1++) 

314 for(x=0;x<nx;x++) 

315 for(y=0;y<nx;y++){ 

316 u = (sqrt(2)/wf)*(x - nx/2 +0.5) * dx; 

317 v = (sqrt(2)/wf)*(y - nx/2 +0.5) * dx; 

318 if (k==0) Hx=l.0; 

319 if (k==l) Hx=2*u; 

320 if (k==2) Hx=4*pow(u,2)-2; 

321 if (k==3) Hx=8*pow(u,3)-12*u; 

322 if (k==4) Hx=16*pow(u,4)-48*pow(u,2)+12; 

323 if (k==5) Hx=32*pow(u,5)-160*pow(u,3)+120*u; 

324 if (k==6) Hx=64*pow(u,6)-480*pow(u,4)+720*pow(u,2)-120; 

325 if (k==7) Hx=128*pow(u,7)-1344*pow(u,5)+3360*pow(u,3)-1680*u; 
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326 if (k==8) Hx=256*pow(u,8)-3584*pow(u,6)+13440*pow(u,4)-13440* 

327 pow(u,2)+1680; 

328 if (k==9) Hx=512*pow(u,9)-9216*pow(u,7)+48384*pow(u,5)-80640* 

329 pow(u,3)+30240*u; 

330 if (1==0) Hy=l.0; 

331 if (1==1) Hy=2*v; 

332 if (1==2) Hy=4*pow(v,2)—2; 

333 if (1==3) Hy=8*pow(v,3)-12*v; 

334 if (1==4) Hy=16*pow(v,4)-48*pow(v,2)+12; 

335 if (1==5) Hy=32*pow(v,5)-160*pow(v,3)+120*v; 

336 if (1==6) Hy=64*pow(v,6)-480*pow(v,4)+720*pow(v,2)-120; 

337 if (1==7) Hy=128*pow(v,7)-1344*pow(v,5)+3360*pow(v,3)-1680*v; 

338 if (1==8) Hy=256*pow(v,8)-3584*pow(v,6)+13440*pow(v,4)-13440* 

339 pow(v,2)+1680; 

340 if (1==9) Hy=512*pow(v,9)-9216*pow(v,7)+48384*pow(v,5)-80640* 

341 pow(v,3)+30240*v; 

342 

343 rsq = ((x-nx/2+0.5)*(x-nx/2+0.5)+(y-nx/2+0.5)*(y-nx/2+0.5))*dx*dx; 

344 phi = ( rsq*(tau-tauw)/(wo*wo*wf*wf)-(k+1+1)*atan((tau-tauw)/zo) ); 

345 

346 Abfinalr[k][1][x][y] = Ao * wo * Hx * Hy * exp(-rsq/(wf*wf)) * 

347 cos(phi)/wf*sqrt(1/(pow(2,k)*pow(2,1)*factorial(k)*factorial(1))); 

348 

349 Abfinali[k][1][x][y] = Ao * wo * Hx * Hy * exp(-rsq/(wf*wf)) * 

350 sin(phi)/wf*sqrt(1/(pow(2,k)*pow(2,1)*factorial(k)*factorial(1))); 

351 

352 Powerfinal[k][1] = Powerfinal[k][1] + 

353 Abf inalr [k] [1] [x] [y] *Abf inalr [k] [1] [x] [y] + 

354 Abf inali [k] [1] [x] [y] * Abf inali [k] [1] [x] [y] ; } 

355 

356 printf("P_00(1) = %g\n",Powerfinal[0][0]); 

357 

358 /* Finding total power in At at tau=l *(should equal the tau=0 case)*/ 

359 for(x=0;x<nx;x++) 

360 for(y=0;y<nx;y++) 

361 PI = PI + Atr[x] [y]*Atr[x] [y] + Ati [x] [y] *Ati[x] [y] ; 

362 

363 printf("P_T(1) = °/„g\n", PI); 

364 

365 /* orthogonality calculation */ 

366 K = (2*dx*dx)/(pi*zo); 

367 for(k=0;k<10;k++) 
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368 

369 

370 

371 

372 

373 

374 

375 

376 


for(1=0;1<10;1++) 
for(x=0;x<nx;x++) 
for(y=0;y<nx;y++) { 
cfinalr[k] [1] = cfinalr[k] [1] + K * 

(Atr [x] [y] *Abf inalr [k] [1] [x] [y] + Ati [x] [y] *Abf inali [k] [1] [x] [y]); 


cfinali [k][1] = cfinali [k] [1] + K * 

(Ati [x] [y] *Abf inalr [k] [1] [x] [y] - Atr [x] [y] *Abf inali [k] [1] [x] [y]); > 


377 /* calculating cmn squared */ 

378 for(k=0;k<10;k++) 

379 for(1=0;1<10;1++) { 

380 c2magfinal[k] [1] = (cfinalr[k][1]*cfinalr [k] [1] + 

381 cfinali[k][1]*cfinali[k][1]); 

382 c2sumFinal = c2sumFinal + c2magfinal[k][1];> 

383 

384 printf("\nNormalized initial cmn amplitudes squared\n"); 

385 for(k=0;k<10;k++) 

386 for(1=0;1<10;1++) { 

387 fS[k][l] = c2magf inal [k] [1]/c2sumFinal; 

388 if (fS[k][l] > 0.01) printf ("Cmn2[%d] [°/ 0 d] = “/„g\n" ,k,l,fS[k] [1]) ;} 

389 printf("\n"); 

390 printf("\nClosing files...\n"); 

391 fclose(input); 

392 fclose(initial); 

393 fclose(final); 

394 fclose(ytau); 

395 


396 } /*end main() */ 


397 

398 int factorial(int a) 

399 { 

400 if ((a==l) II (a==0)) 

401 return 1; 

402 else { 

403 a *= factorial(a-1); 

404 return a; } } 

405 /* End of Cmn.c */ 
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2. CMN.IN - INPUT FILE 

/* Begin input file */ 


1 

1 0 

2 


2 

1 2 

0 


3 

0 0 

2 


4 

0 4 

0 


5 

0 0 

4 


6 

0 2 

7 


7 

0 4 

9 


8 

0 3 

8 


9 

0 5 

3 


10 

0 0 

6 


11 

1.0 

10 

100 

12 

0.5 

l.i 

0 400E-6 

13 

*/.c0 

mO 

nO 

14 

*/.cl 

ml 

nl 

15 

*/.c2 

m2 

n2 

16 

°/oC3 

m3 

n3 

17 

°/oC4 

m4 

n4 

18 

°/.c5 

m5 

n5 

19 

%c6 

m6 

n6 

20 

y.c7 

m7 

n7 

21 

%c8 

m8 

n8 

22 

°/«c9 

m9 

n9 

23 

°/oAo 

W nx 

24 

°/oZO 

tauw dtau 

/ * End 

input : 

file */ 


3. COEFFS.M - MATLAB SCRIPT 

1 close all 

2 clear 

3 clc 

4 load initial.out 

5 load final.out 

6 load Cmn.in 

7 load ytau.out; 

8 
9 

10 ytau = ytau’; 

11 initial = initial’; 

12 final = final’; 
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13 

14 

15 W = Cmn(ll,2); 

16 nx = Cmn(ll,3); 

17 dtau = Cmn(12,3); 

18 tau=0:dtau:1-dtau; 

19 x = 0:1:nx-l; 

20 

21 [npx,timesteps] = size(ytau); 

22 dtau = 1/timesteps; 

23 tau = 0:dtau:1-dtau; 

24 x = 0:l:npx-l; 

25 

26 

27 [x_initial,y_initial] = meshgrid(0:W/nx:W-W/nx); 

28 [x_final,y_final]=meshgrid(0:W/nx:W-W/nx); 

29 

30 

31 figure(’Color’,’white’) 

32 set(gca,’FontSize’,25) 

33 pcolor(x_initial,y_initial,initial) 

34 colorbar 

35 axis square 

36 shading flat 

37 shading interp 

38 title(’initial waveform’) 

39 xlabel(’x’) 

40 ylabel(’y’) 

41 colormap(hot) 

42 

43 figure(’Color’,’white’) 

44 set(gca,’FontSize’,25) 

45 surf(initial) 

46 title(’initial surface plot’) 

47 xlabel(’x’) 

48 ylabel(’y’) 

49 zlabel(’a_i*a_i’) 

50 colormap(hot) 

51 

52 figure(’Color’,’white’) 

53 set(gca,’FontSize’,25) 

54 pcolor(x_final,y_final,final) 
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55 xlabel(’x’) 

56 ylabel(’y’) 

57 colorbar 

58 axis square 

59 shading flat 

60 shading interp 

61 title(’final waveform’) 

62 colormap(hot) 

63 

64 figure(’Color’, ; white’) 

65 set(gca,’FontSize’,25) 

66 surf(final) 

67 xlabel(’x’) 

68 ylabel(’y’) 

69 zlabel(’a_f*a_f*’) 

70 title(’final surface plot’) 

71 colormap(hot) 

72 

73 figure(’Color’,’white’) 

74 set(gca,’FontSize’,25) 

75 [Tau,X]=meshgrid(tau,x); 

76 pcolor(Tau,X,ytau) 

77 shading flat 

78 shading interp 

79 title(’Laser amplitude y vs. 

80 xlabel(’\tau’) 

81 ylabel(’x’) 

82 colormap(hot) 


\tau’) 
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APPENDIX B. USEFUL RELATIONS 


Quick reference for some important equations 


1. FREE ELECTRON LASER OVERVIEW 


Resonance Condition 


Gain 


Extraction 


A = 

G = 

V = 



P-Po 

Pa 


_ Optical power _ 

Initial electron beam power 


Equation (II. 1) 
Equation (II.2) 
Equation (II.3) 


2. FREE ELECTRON LASER THEORY 


Pendulum Equation 


Power 


FEL wave equation 


C= M cos(c + </>) 

P=\a \ 2 



d_ 

dr 


a(r,r) 


=< — je ^ > 


Equation (III.32) 
Equation (III.37) 
Equation (III.64) 
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3. OPTICAL THEORY 


F. Trial Solution 

-G+TL) 
a = a 0 e v 

Equation (IV.3) 

F. Solution 

r 2 % ( \ T —arctan T- 

a = a 0 —e ^eVv 2 

° W 

) Equation (IV. 20) 

HOM Trial Solution 

( f2 'i 

a = a 0 gheP (p+I ^' 

Equation (IV. 24) 

HG Mode Solution 

a m , n = a 0 ^H m H n eT^^ 

Equation (IV.45) 

HOM Phase 

<j>m,n = J o 2 ^ 2 (m + n + 1) arctan I Equation(IV.46) 

Normalized HG mode 

A — 1 

^m,n s '2<"2" ,n\n' 

Equation (IV.60) 

Total field HG 

At = Ylm,n Cm,nA m , n 

Equation (IV.61) 

Normalized Coefficients 

Cm.n = na 2 w 2 f :0C / ^ -4t -4 m . n dxdy Equation (IV.66) 

SIMULATIONS 



Total Power in optical field Pt = § a 2 0 w 2 0 Ym n ^'m.n 

Equation (V.13) 

Fractional Power Spectrum S[m][n] = 

Equation (V.15) 
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